intro

Transcriptomics

Transcriptomics is de analyse van de RNA transcriptie’s geproduceerd door het genotypen op een bepaald moment wat een link geeft tussen het genoom, het proteoom en het cellulaire phentype. Het is een stap in de Omics pipeline waar dingen zoals ziekten en genetische afwijkingen en het effect er van word onderzocht. Transcriptomics in ons onderzoek word gedaan door de RNA transcripties te vergelijken van 2 groepen muizen. Een jonge groep muizen van 8-10 weken en een oudere groep muizen van 100-120. In ons onderzoek kijken naar de expressie in een paar gen pathways: 04662, 04660, 04213, 04211, 04640.

De 04662 is de B cell receptor signaling pathway deze laat zien ““.

De 04660 is de T cell receptor signaling pathway deze laat zien ““.

De 04211 is de Longevity regulating pathway deze laat zien ““.

De 04640 is de Hematopoietic cell lineage deze laat zien ““.

In de gene pathways kan er worden gekeken naar de verandering in expressie, hier is dus te zien of de gevonden variant mutatie effect heeft op de expressie van het gen. Een gen pathway ziet er uit als volgt.

Figuur 1: Voorbeeld gene pathway
Figuur 1: Voorbeeld gene pathway

Hier is een door ons met KEGG + pathview library gemaakte gene pathway te zien. Rechts boven is een kleuren schaal te zien met -1, 0 en 1. De kleuren schaal en cijfers betekenen ook wat, wanneer de gen een witte/grijze kleur heeft is de verandering in de expressie 0 wat dus betekent dat de expressie van dit gen normaal is. Wanneer een gen block een groene kleur heeft en dus richting de -1 gaat zegt dit dat de expressie is minder is dan de norm. Wanneer een gen een rode kleur heeft en dus richting de 1 gaat zegt dit dat het een overexpressie heeft meer dan het normaal.

De verschillen in deze expressies hebben veel invloed op de werking van de pathway.

Hier is een flowchart te zien die laat zien hoe wij dit project gaan doen. We beginnen dus met de data downloaden, controleren de kwaliteit en trimmen. Wanneer de kwaliteit goed is gaan we readmappen met HISAT2 daarna data met samtools sorteren en indexeren. De Data wordt dan verwerkt met Feature Counts om de gen expressie vast te stellen. Dan word de data geanalyseerd op 2 manieren. Met de wilcoxon rank sum test, deze test beoordeelt niet alleen het teken van de verschillen, maar bepaalt ook de grootte van de waargenomen verschillen in de gegevens. De andere test is DEseq2, Deseq gebruikt een statistisch model om uitterekenen wat de verschillen in gen expressie zijn tussen 2 of meer groupen samples (RNA). Aan de hand van de feature count data en de DEseq2 uitkomt worden visualisaties gemaakt.

KEGG - is een knowledge base voor systematische anlyse van gen functies, waar genomics informatie gelinked kan worden met functionele info. Er zijn ook visualisatie die hier uit kunnen komen zoals pathways waar over en onder expressie van genen word laten zien in een gen pathway

Volcano plot - is een plot gemaakt door het plotten van de negatieve logartime van de p value op de y axis meestal met een base van 10. Met de log fold change op x axis. Een volcano plot laat de statistische significantie zien met P value tegen de groote van de fold change.

figuur 2: Flowchart - transcriptomics
figuur 2: Flowchart - transcriptomics

Hier zijn alle library’s die wij gebruiken in dit project.

## Loading required package: S4Vectors
## Loading required package: stats4
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
##     tapply, union, unique, unsplit, which.max, which.min
## 
## Attaching package: 'S4Vectors'
## The following object is masked from 'package:utils':
## 
##     findMatches
## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname
## Loading required package: IRanges
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## Loading required package: ggplot2
## Loading required package: ggrepel
## 
## Attaching package: 'data.table'
## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift
## The following object is masked from 'package:GenomicRanges':
## 
##     shift
## The following object is masked from 'package:IRanges':
## 
##     shift
## The following objects are masked from 'package:S4Vectors':
## 
##     first, second
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:S4Vectors':
## 
##     expand

Code - path melding

Wanneer er geen input path staat bij het BASH script ga er dan vanuit dat de code in de map is uitgevoerd waar de input bestanden staan.*

SRA download

De SRA’s van de RNA seq worden gedownload met het onderstaande stuk code. De gebruikte accention list staat in de github repository.

De eerste stap is het downloaden van de genetische data die verkregen was door het sequencen van de samples. Met behulp van de SRA run selector van NCBI is een selectie van SRA’s gemaakt die gedownload moeten worden.

Een SRA (Sequence Read Archive) is een gecomprimeerd archief dat de sequencing reads bevat.

Om deze bestanden te downloaden wordt gebruikt gemaakt van “prefetch”, deze commandline tool is onderdeel van de SRA toolkit.

In de volgende stap worden deze bestanden uitgepakt.


prefetch $(</students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt) \
--output-directory "/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/" --max-size 200G

Gebruikte samples:

In ons onderzoek gebruiken we de 2 muizen groepen, de oude en jonge muizen. De oude muizen zijn 100-120 weken oud en de jonge muizen zijn 8-10 weken oud. Deze twee groepen hebben allebei vier replicaten. De samples van alle muizen zijn verkregen met een Illumina NextSeq 500.

Oude muizen: (Genotype: Rag2-/-)

SRR21754423
SRR21754422
SRR21754421
SRR21754420

Jonge muizen: (Genotype: Rag2-/-)
SRR21754408
SRR21754417
SRR21754418
SRR21754419

Kwaliteits check

De gedownloade bestanden zitten nog in hun SRA format en dit is niet de filetype die wij gebruiken voor onze tools.

Als eerste moet de kwaliteit van de data gecheckt worden en dit wordt gedaan met het programma FastQC dit programma leest de ASCII characters uit de file die bij elke sequentie staan waardoor het weet wat de kwaliteit is van de data. De voledige werking van FastQC word later nog wat dieper uitgelegt.

Om de bestanden van .SRA om te zetten naar .fastq word fasterq-dump gebruikt. Deze kijkt in de SRA files en pakt ze uit naar hun fastq format. Fasterq-dump is deel van de SRA toolkit (onze gebruikte versie is 3.1.1)

Een FastQ file ziet er als volgt uit.

De file heeft een title waar de sample naam word weergegeven + de lengte. Daarna is de sequentie te zien, Dan is de titel te zien van de bijbehoorende kwaliteit sequentie. De laatste regel is dus de kwaliteits score sequenties deze is in ASCII

De kwaliteit word dus gescoord in ASCII characters, in het figuur hier onder is de score te zien van een voorbeeld sequentie, ook is te zien dat de score van het ! teken de laagst mogelijk kwaliteit is. De K is de hoogste kwaliteit mogelijk.

De formule die hier wordt gebruikt is de Phred kwaliteits score, deze score is logaritmisch gerelateerd aan de waarschijnlijkheid dat de “base call” verkeerd is.

Q = -log(E) waarbij Q de phred score is en E de waarschijnlijkheid van verkeerde base call.

De illumina gesequencenste sequenties hebben een afwijking van 33 hier word dus altijd -33 gedaan voordat de score wordt berekent.

Figuur 3: fasta file voorbeeld Bron: https://gencoded.com/index.php/2020/05/20/fastq-format-an-overview/

Bron: https://gatk.broadinstitute.org/hc/en-us/articles/360035531872-Phred-scaled-quality-scores

# Cat leest elke File-parameter in volgorde van Acc_list_transcriptomics.txt  en daarmee kijkt fasterq-dump welke files het moet uitpakken. Met parallel worden meerder files tegelijk verwerkt

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | \
  parallel fasterq-dump -O /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/FASTQ {}
  # hierboven wordt de map output neergezet

Kwaliteitscontrole met fastqc

Voor het trimmen: FastQC wordt gebruikt voor het trimmen. Dit programma kan de scores lezen van de fastQ files en geeft hier een visuele output van. Dit zijn verschillende grafieken en tabellen hier onder is een tabel met uitleg voor elke visualisatie / tabel.

Catagorie Uitslag Hoe lees je dit
Basic s tatistics Hier krijg je een tabel met verschillende dingen zoals de Filename, file type, sequence length, %GC, Total sequences, Encoding en de sequences flagged as poor quality Het aflezen van deze tabel is lezen en checken of er geen fouten in zitten.
Per base sequence quality Een grafiek zoals te zien is in figuur 1

De grafiek heeft 3 vlakken die van verschillende grote en kleuren. De kwaliteit ranged van 0-40, 28-40 is goed and gekleurd groen, 20-28 is gekleurd geel/oranje, 0-20 gekleurd rood.

Verder representeert de gele box de 25th to 75th percentile. De zwarte lijnen geven de 10th en 90th percentile weer. De blauwe lijn geeft de gemiddelde scores voor kwaliteits controle score voor de nucleotide. gebaseerd op deze dingen is te zien dat de base van 1 tot 100 goede kwaliteit hierna gaat de kwaliteit sterk naar beneden. wat al aan geeft dat er getrimmed moet worden.

Per tile sequence quality Hier is een plot te zien waar de reads worden weergegeven als een soort heatmap. waarbij elke tegel 1 een read weergeeft op een bepaalde positie deze kwaliteit van de reads worden vergeleken met elkaar. Het aflezen van deze grafiek is door te kijken naar de kleur per tegel per positie. een donkerblauwe tegel betekent goede kwaliteit en hoe lichter de tegel word hoe slechter de kwaliteit dus lichtblauw betekent een slechte kwaliteit. Op de Y-as staat dan in welke tile het is en de X-as welke positie.
Per sequence quality score Een grafiek die je de gemiddelde kwaliteit op de
x-as en de nummer van sequences met gemiddelde op de y-as.
Het beste is wanneer de meeste reads een hoge gemiddelde kwaliteits score hebben en er geen grote dip in de grafiek is, dit betekent een lage kwaliteit.
Per base sequence content Hier is een plot te zien met een y-as waar 0-100 aangegeven wordt en een x-as met de “positie in read (bp)” in de grafiek staan 4 lijnen met het percentage per base Het volgende is uit de grafiek te halen er is significante variatie in de nucleotide distributie aan het bein van de reads positie 1-10. Dit zou kunnen zijn door de voorbereiding of de vooroordelen in het sequencen. A, T, C en G zijn niet gelijk gerepresenteerd. na de 10de positioe zijn de base wat gestabilizeerd wat aan geeft dat de sequence kwaliteit in de rest van de reads hoger zijn. Het afgekeurde kruis komt dus door de choas van 1-10.
Per sequence GC content Per sequence GC content geeft weer een plot weer met 2 lijnen. Een blauwe lijn die de theoretische distrubutie aangeeft wat dus een richtlijn is, en de GC count per read wat dus de gelezen data is. Het aflezen wordt door de twee lijnen vergelijken. het is ideaal als de twee lijnen overlappen of dichtbij elkaar liggen. Wanneer er meerdere pieken zijn die afwijken van de blauwe theoretische lijn kan dat betekenen dat er misschien sprake is van besmetting of sequencing errors. Er komt dus een rood kruis wanneer de GC abnormaal is vergeleken met de theoretische verwachting wat dus zegt da de algemene kwaliteit niet goed is.
Per base N content De per base N content grafiek geeft de frequentie weer van “N” basecalls op elke positie in de reads. De “N” staat voor een onzekere of onbekende nucleotide, deze kon de sequencer niet identificeren als 1 van de base (A, T, C, G) X-as geeft de positie weer in de reads, Y-as geeft het percentage van reads met een “N” base op elke positie. Een hoge waarde betekent dat de sequencer op die positie vaak onzeker was over welke nucleotide er aanwezig was. De verwachting is dat er een zeer laag percentage N’s in de sequence zit, de standaard hiervoor is <1%. afwijkende resultaten kan wijzen op problemen zoals slechte kwaliteit van de sequentie. Ze komen vaak voor aan het begin of eind van de reads.
Sequence length dis tribution De sequence length distribution grafiek laat de verdeling van de lengtes van de reads zien. X-as geeft de lengte van de sequenties (in basenparen) en de Y-as toont aantal reads van die specifieke lengte. Deze grafiek is belangrijk bij NGS omdat afwijkingen in de sequentie lengte kunnen wijzen op fragmentatie of technische fouten tijdens sequencering . Een ideaal resultaat is een scherpe piek op 1 specifieke lengte wanneer je 150 bp doelreadlengte hebt zou de meeste sequencing output op precies 150 bp moeten vallen wat een piek bij die lengte zou moeten opleveren. Als er meerdere pieken zijn of een vrede spreiding van lengtes kan dit beteken dat er sequencing fouten, slechte adaptertrimming of degradatie van het DNA-monster.
Sequence du plication levels

De Sequence duplication levels grafiek laat het percentage van sequences zien die meer dan 1x voorkomen. Duplicaties kunnen voorkomen door technische “artifacts” tijdens de sequencing en andere factoreren zoals P CR - amplificatie, en kunnen de diversititeit van de dataset verminderen. Dit heeft invloed op downstream analyses zoals genexpressie of va r iantdetectie.

X-as laat het aantal keren zien dat de specifieke sequentie wordt gedupliceerd

Y-as laat het percentage van reads dat voorkomt met duplicatie niveau zien

Het aflezen van de grafiek is door te kijken naar de twee lijnen in de grafiek de Duplicated sequences (meestal een rode lijn) geeft het percentage met unieke sequenties weer zonder correctie voor natuurlijke duplicaties. Het geeft een beeld van hoeveel van de reads in de dataset meerdere keren voorkomen zonder te differentiëren tussen technische en biologische duplicates. in de grafiek wil je graag een scherpe daling zien waarbij de meeste reads een duplicatie niveau van 1 hebben en het percentage gedupliceerede sequenties daarna snel afneemt.
De total sequences lijn (meestal een blauwe lijn) corrigeert voor verwachte natuurlijke duplicates deze laat zien hoe de duplicatie eruit zou zien zonder technische artefacten en biedt een eerlijker beeld van hoeveel sequenties overgedupliceerd zijn de de sequencering zelf.
Als er een groot verschil is tussen de lijnen betekent dit dat er duplicatie is onstaan door technische factoren zoals PCR duplicatie in plaats van biologische oorzaak. Te hoge aantal dupolactie kan probleem zijn voor downstream analyses en je wil dus dat de lijnen dichtbij elkaar liggen.
Overre presented sequences Dit is een tabel die sequenties van op zijn minst 20bp die vaker voorkomen dan 0.1% van de totale nummer van sequenties. In de tabel staat de sequenties u itgeschreven, de count, het percentage en de w a arschijnlijke bron Wanneer er een abnormaliteit te zien was in de Per sequence GC content grafiek kan er in deze tabel worden gekeken om de bron te vinden. als het niet staat als een bekende adapter of “vector”, kan het helpen om de data te blasten o te identiteit te vinden in de tabel.
Adapter content Een grafiek waar verschillende adapters staan Het geeft aan of de sequenties adapterfragmenten bevatten en van verschillende apparaten. als deze aanwezig zijn is het te zien door af te lezen in de grafiek welke positie er zijn om ze vervolgens weg te trimmen. deze “adaptercontent” is er voor identificatie van de DNA

Tabel 1: Uitleg fasta uitlezing

Het volgende commando is uitgevoerd in de FASTQ/ directory.

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
fastqc -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/fastqc_output/voor_trimmen/ {}_1.fastq {}_2.fastq 

Vervolgens met met “multiqc .” in de output directory een multiqc rapport gemaakt van de resultaten van fastqc.

MultiQC wordt gebruikt om naar alle fastqc output files te kijken. Multiqc zoekt een directory voor FastQC output en compileert deze files zodat de algemene resultaten over meerdere files beter te begrijpen is. De grafieken worden bijvoorbeeld over elkaar gelegd om de verschillen te zien.

Figuur 4: multiqc sequentie kwaliteit
Figuur 4: multiqc sequentie kwaliteit

Na analyse van de fastQC is te zien dat onze sequence reads korte sequenties zijn van 150 bp. Deze corresponderen met een deel van een DNA sequentie. Het genoom kan niet in zijn geheel gelezen worden door technische beperking van de gebruikte illumina machines. Hoe langer de reads hoe meer kans dat er fouten gemaakt worden tijdens het sequencen.

daarom als er genoeg reads zijn die het laatste stuk bevatten kan er een stuk afgeknipt worden. Ook moeten de adapters verwijderd worden, adapters zijn korte stukken DNA die op de flow cell zitten als linkers. bron: https://www.lubio.ch/blog/ngs-adapters

Voor het trimmen wordt gebruik gemaakt van Trimmomatic ( versie 0.39 ) Trimmomatic gebruikt Maximum Information Quality Filtering om te bepalen of een read moet worden getrimmed + de paramenters van de gebruiker. Maximum Information Quality Filtering probeert de formule te bepalen bij welke informatie de hoogste kwaliteit heeft door verschillende eigenschappen van de informatie te combineren, zoals: Relevantie (afgeleid via de sigmoidfunctie), Correctie of betrouwbaarheid (door het product van correctieprobabiliteiten). Het doel is om de score te maximaliseren voor informatie die zowel relevant als zeer betrouwbaar is.

De door ons opgegeven parameters staan in het codeblock deze zijn.

MINLEN:40 dit zegt dat de minimale lengte van de reads 40 is.

SLIDINGWINDOW 4:20 het eerste nummer specificeert de grootte van de sliding winden en het tweede nummer is de gemiddelde read kwaliteit binnen de window van het eerste nummer. Een sliding window is dat er steeds in dit geval 4 basen bekeken worden en dat keer op keer.

ILLUMINACLIP is het path naar een bestand dat adapter sequenties bevat voor illumina adapters. De gebruikte adapter: TruSeq3-Se.fa:2:30:10 De 2:30:10 betekenen het volgende: 2 is de “seed mismatch”, het aantal mismatches dat toegestaan is in een sequentie die een adapter kan zijn.30 is de “palindrome clip threshold en 10 is de”simple clip threshold”, specificeert hoe accuraat de match tussen de adapter sequentie en de mogelijke adapter in de read.

Hier is een foto te zien van hoe Trimmomatic weet wanneer hij moet knippen How Maximum Information mode combines uniqueness, coverage and error rate to determine the optimal trimming point

Trimmomatic neemt 3 scores, de length Threshold, coverage factor en de Error rate factor en vormt hier 1 totale score is en trimt hem op de piek zoals te zien in de grafiek hierboven.

Dit is gebaseerd op deze formule.

 formula

Deze formule berekent een score \(\text{Score}(l)\) door een combinatie van drie elementen: een sigmoid-achtige factor om de score op basis van \(l\) te dempen, een machtsverhouding van \(l\) met een parameter \(s\), en het product van correctiewaarden \(P_{\text{corr}}[i]\), die allemaal samen de uiteindelijke score beïnvloeden.

Bron: http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'TrimmomaticPE -threads 80 ' \
    '/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/FASTQ/{}_1.fastq' \
    '/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/FASTQ/{}_2.fastq' \
    '/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_forward_paired.fastq' \
    '/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/unpaired/{}_forward_unpaired.fastq' \
    '/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_rev_paired.fastq' \
    '/students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/unpaired/{}_rev_unpaired.fastq' \
    'ILLUMINACLIP:/students/2024-2025/Thema05/3dconformatieChromatine/Trimmomatic/Trimmomatic-0.39/adapters/TruSeq3-SE.fa:2:30:10' \
    'MINLEN:40 ' \
    'SLIDINGWINDOW:4:20'

Kwaliteitscontrole na trimming met fastQC

Na het trimmen van de reads kijken we weer naar de kwaliteit dit is om te kijken of we de parameters moeten aan passen van de trimmomatic of dat we met deze getrimmde data kunnen werken.

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
fastqc -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/fastqc_output/na_trimmen/ {}_forward_paired.fastq {}_rev_paired.fastq 

Indexeren referentie genoom

Om te kunnen zien of onze reads verschillen met de normale expressie is het belangrijk dat wij een referentie genoom meenemen in onze analyse. De tools die wij gebruiken hebben tevens het referentie genoom nodig zodat ze functioneren.

Door middel van readmapping worden de reads dus vergeleken met het bekende genoom, in ons geval de mm39 het op dit moment nieuwste muis referentiegeoom. De reads worden er dus tegen aan gelegd zodat hun locatie in het genoom bekend wordt.

Het programma dat gebruikt gaat worden voor het alignen is HIsat2, deze heeft dus zoals eerder al gezegd een referentie genoom nodig maar deze moet dus geindexeerd zijn en deze moet eerst gemaakt worden.

Dit kan gedaan worden met Hisat2 build. Het referentie genoom mm39 word gebruikt ook wel met de bestand naam GCF_000001635.27, In het bestand met het referentie genoom zijn twee referentie genomen aanwezig GCF_000001635.27 en GCA_000001635.9/. Het verschil is dat GCF van Refseq is en GCA van GeneBank.https://www.ebi.ac.uk/training/online/courses/functional-genomics-ii-common-technologies-and-data-analysis-methods/rna-sequencing/performing-a-rna-seq-experiment/data-analysis/read-mapping-or-alignment/

Bron: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3375638/
Bron Dowload ref genoom : https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000001635.27/

In de onderstaande code chunk wordt het geïndexeerde muisgenoom “hisat2_index” genoemd.Daarna worden twee samples (forward read en reverse read) ge-aligned met het zojuist verkregen geïndexeerde referentiegenoom.Deze twee bewerkingen worden beiden met HIsat2 uitgevoerd.

HISAT2 is een gevoelig en snel alignment programma voor het in kaart brengen van ngs reads dit is voor zowel DNA als RNA naar een populatie van menselijke genomen en naar een enkel referentiegenoom. Gebaseerd op een uitbreiding van BWT voor grafieken (Sirén et al. 2014) is er een graph FM index (GFM) ontworpen en geïnplementeerd, een originele aanpak en de eerste implementatie ervan. Naast het gebruik van één globale GFM index die een populatie van menselijke genomen vertegenwoordigt, gebruikt HISAT2 een grote set van kleine GFM indexen die gezamelijk het gehele genoom bestrijkt. Deze kleine indexen ook wel lokale indexen genoemd, gecombindeerd met verschillende uitlijningsstragieën, maken snelle en nauwkeurige alignments van sequencing reads mogelijk. Dit nieuwe indexeringsschema wordt een Hierarchical Graph FM index (HGFM) genoemd.

#Hisat build wordt dus gebruikt om de ref genoom correct te aligneren de -p optie laat HISAT2 een specifiek aantal parallel zoek threads gebruiken die 60 zijn door ons aangegeven
hisat2-build -p 60 /students/2024-2025/Thema05/3dconformatieChromatine/Mapping_ref/ncbi_dataset/ncbi_dataset/data/GCF_000001635.27/GCF_000001635.27_GRCm39_genomic.fna hisat2_index

Read mappig:

Hisat index

Hier wordt met hisat2 de sequenties ge aligned aan het zonet gemaakte referentie genoom. De -p optie laat HISAT2 een specifiek aantal parallel zoek threads gebruiken die 80 zijn door ons aangegeven. -x geeft de basename van de index voor het reference genoom, in ons geval de hisat2_index. -1 geeft door komma’s gescheiden lijst met bestanden die mate 1’s bevatten (bestandsnaam bevat meestal _1), bijvoorbeeld -1 flyA_1.fq,flyB_1.fq. Sequenties die met deze optie worden opgegeven, moeten bestand-voor-bestand en lees-voor-lees overeenkomen met die welke zijn opgegeven in <m2>. Lezen kan een mix van verschillende lengtes zijn. Als - is opgegeven, leest hisat2 de mate 1’s van de “standaard in” of “stdin”-bestandshandle.

Daarna leest hisat de forward- en reverse paired reads de {} geeft aan dat de bestand naam alles mag zijn maar het moet _forward_paired.fastq of _rev_paired.fastq bevatten om een file te zijn die mee word genomen in de hisat.

-S zegt dat Bestand om SAM-uitlijningen naar te schrijven. Standaard worden uitlijningen naar de “standard out” of “stdout” filehandle geschreven (d.w.z. de console).

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
hisat2 -p 80 -x hisat2_index -1 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_forward_paired.fastq -2 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/trimmomatic_output/paired/{}_rev_paired.fastq -S /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/hisat2_mapped/{}.sam &

De output van de read mapping is een sam bestand, een sam bestand ziet er zo uit:

Figuur 7: structuur sam bestand
Figuur 7: structuur sam bestand

Hier is het format te zien van een sam file. Er staan verschillende dingen in zoals de alignement info, de sequentie read, mapping quality read en nog verschillende dingen die belangrijk zijn voor de tools in de pipeline.

Sorteren en indexeren van .sam

Omzetten van .sam naar .bam: De reden om de file van een .sam naar een .bam om te zetten is omdat .bam files beter manipuleer baar zijn en ook sneller dus daarom beter zijn om mee te werken.

#Cat leest elke File-parameter in volgorde van 
cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
samtools view -@40 -b /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/hisat2_mapped/{}.sam > /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/{}.bam &

Hier gebruiken we samtools view dit zet de .sam file over naar een .bam, de -@ optie word gebruikt om meerdere threads te selecteren in ons geval 40 de -b optie geeft aan welke soort file wij als output willen en -b staat natuurlijk gelijk aan .bam De andere twee lines geven de path aan waar de input files staan, en waar de output heen moet. & geeft aan dat het programma in de achtergrond moet runnen.

Index maken: Indexeer coördinaten-gesorteerde BGZIP-gecomprimeerde SAM-, BAM- of CRAM-bestanden voor snelle willekeurige toegang.Deze index is nodig wanneer regio-argumenten worden gebruikt om de samtools-weergave en soortgelijke opdrachten te beperken tot bepaalde regio’s van belang.

samtools index roept natuurlijk de index aan, de -@ optie word gebruikt om meerdere threads te selecteren in ons geval 80 en alles wordt in parallel gerund dus meerdere opdrachten tegelijk.

Verder is nog het output path te zien, hier worde

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
samtools index -@80 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/sorted_{}.bam
  

Samtools sort:

Sorteer uitlijningen op meest linkse coördinaten, op gelezen naam wanneer -n of -N worden gebruikt. De gesorteerde uitvoer wordt standaard naar de standaarduitvoer geschreven of naar het opgegeven bestand (out.bam) wanneer -o wordt gebruikt.

samtools sort roept dus de sorting aan. @ zorgt dat er een bepaald aantal threads worden gebruikt in ons geval 80. -o geeft aan dat het ouput bestand een bepaald bestand type moet zijn dit is bij ons BAM . -n geeft aan dat ons bestand op naam gelezen word.

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
samtools sort -@80 -O BAM -n  /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/{}.bam -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/sorted_{}.bam

samtools fixmate

samtools fixmate – vult in mate coördinaten en insert “size fields” -m voegt mate scored toe en deze zijn nodig bij markdup om de beste reads te behouden -@ 80 zegt dat er 80 threads moeten worden gebruikt voor deze opdracht verder zijn er nog 2 paths die de input en output geven.

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'samtools fixmate -m -@80 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/sorted_{}.bam /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_{}.bam'

Sorteren op coördinaten, dit is de standaard dus hoeft er geen argument voor meegegeven worden.Vraag: Waarom zijn de gesorteerde bestanden kleiner dan de niet gesorteerde?

Sorteer uitlijningen op meest linkse coördinaten, op gelezen naam wanneer -n of -N worden gebruikt. De gesorteerde uitvoer wordt standaard naar de standaarduitvoer geschreven of naar het opgegeven bestand (out.bam) wanneer -o wordt gebruikt.

samtools sort roept dus de sorting aan. @ zorgt dat er een bepaald aantal threads worden gebruikt in ons geval 80. -o geeft aan dat het ouput bestand een bepaald bestand type moet zijn dit is bij ons BAM .

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'samtools sort -@80 /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_{}.bam -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_coords_{}.bam'

Duplicaten verwijderen: samtools markdup markeert de duplicaten in een gecoördineerd en gesorteerde file

samtools markdup roept de markdup functie aan in samtools. -@ 80 zegt dat er 80 threads moeten worden gebruikt voor deze opdracht. -r zegt dat de gemarkeerde duplicates moeten worden verwijderd. -s zegt dat er wat basis statistieken er bij moeten worden gevoegd.

cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'samtools markdup -@80 -r -s /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/fixed_mates_sorted_coords_{}.bam /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/dedup_{}.bam'

Exploratory Data Analysis:

Featurecounts:bron: https://subread.sourceforge.net/featureCounts.html-T 64 was het maximale aantal threads dat het programma kan gebruiken.

featureCounts roept het programma feature count aan.

-T geeft de aantal threads aan die het programma mag gebruiken, in dit geval 64 dit is de maximale aantal threads dat het programma kan gebruiken door het limiet van feature counts. dit is verder te lezen in de bron hier boven.

-B creërt 1 blok index. De bouw index zal niet worden gesplitst in meerdere stukken. Hoe meer blocken een index heeft, deste langzamer de mapping snelheid. Deze optie zal de -M overschrijven als deze ook word meegegeven.

-p specifeert dat de data paired-end reads zijn, feature counts zal stoppen wanneer het een single-end reads ook wel anders dan gespecificeerd. Om de fragmenten te tellen voor paired -end reads, zal de —countReadPairs parameter moeten worden gespecificeerd

-t dit specifeert de type input van sequencing data. er zijn verschillende value’s die de volgende bevatten, 0 = RNA-seq, 1 = DNA-seq data. De gebruiker moet een bepaalde waar specificeren dit kan echter ook een karakter value zijn. Bij ons is exon aangegeven.

-g specificeert de attribute type gebruikt voor het grouperen van features (eg. exons) naar meta-features (eg. genes) wanneer GTF annotatie “provided”. gene id is de default. Dit attribute type is meestal de gen identifier. Dit argument is handig voor meta feature level samenvatting

-a geeft aan dat er een gen annotatie file is die de chorosomale coördinaten van de exonen van elk gen bevat. De path hiervan komt dan ook na de -a.

-o specifeert de base naam van de index die gemaakt word.

verder zijn dus de paths te zien.


cat /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/SRA/Acc_list_transcriptomics.txt | parallel \
'featureCounts -T 64 -p --countReadPairs -B -t exon -g gene_id -a /students/2024-2025/Thema05/3dconformatieChromatine/snpEff/snpEff/data/mm39/genes.gtf -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/featurecounts_output/{}_counts.txt /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/dedup_{}.bam'

# Test:
featureCounts -T 64 -p --countReadPairs -a /students/2024-2025/Thema05/3dconformatieChromatine/snpEff/snpEff/data/mm39/genes.gtf -o /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/featurecounts_output/count.out /students/2024-2025/Thema05/3dconformatieChromatine/transcriptomics/samtools_output/dedup_*

Resultaten van FeatureCounts:

A s s i g n e d 7 , 4 7 8 , 2 7 5 1 1 , 5 1 9 , 0 3 3 4 0 1 , 9 1 9 1 , 7 9 8 , 7 8 3 1 , 0 3 4 , 2 0 6 1 , 8 1 5 , 6 8 5 2 , 1 2 3 , 1 1 7 1 , 5 5 9 , 9 6 5
U n a s s i g n e d

_ U n m a p p e d
1 0 1 , 5 6 9 9 1 , 9 5 7 3 0 3 , 1 4 7 5 0 1 , 9 1 8 1 , 9 0 0 , 7 6 5 2 , 0 0 6 , 3 6 6 2 , 3 0 5 , 5 5 2 3 , 1 5 3 , 7 0 0
U n a s s i g n e d
R e a d
T y p e
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ S i n g l e t o n
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ M a p p i n g Q u a l i t y
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ C h i m e r a
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ F r a g m e n t L e n g t h
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ D u p l i c a t e
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ M u l t i M a p p i n g
1 , 2 7 8 , 1 8 4 1 , 2 7 8 , 6 6 7 1 , 4 2 3 , 5 8 5 2 , 4 4 4 , 3 4 2 1 , 7 2 8 , 7 3 8 1 , 7 1 8 , 8 6 9 1 , 8 6 8 , 8 4 9 1 , 9 7 3 , 5 7 6
U n a s s i g n e d

_ S e c o n d a r y
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ N o n S p l i t
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ N o F e a t u r e s
7 2 3 , 9 8 2 9 9 9 , 0 0 1 5 8 6 , 7 1 9 1 , 2 8 4 , 5 7 6 4 2 4 , 9 2 5 5 7 4 , 5 0 4 5 6 9 , 8 1 1 7 4 4 , 2 9 8
U n a s s i g n e d
O v e r l a p p i n g
L e n g t h
0 0 0 0 0 0 0 0
U n a s s i g n e d

_ A m b i g u i t y
3 0 , 5 7 5 3 2 , 5 9 8 3 1 , 9 2 3 1 7 , 5 9 8 1 2 , 4 0 6 1 5 , 2 7 6 1 8 , 9 3 9 1 4 , 7 8 3

Table 2: feature count resultaten

Assigned (A): Dit zijn de reads die correct zijn toegewezen aan genen. Hoge waarden hier wijzen op een succesvolle toewijzing en een efficiënte sequencing-run.

Unassigned: Er zijn meerdere subcategorieën hier die de oorzaken van de niet-toewijzing weergeven:

  • Mapping Quality: Reads die door lage mapping-kwaliteit niet konden worden toegewezen.

  • Ambiguity: Reads die aan meerdere genen kunnen worden toegewezen en daardoor ambigu zijn.

  • Multi Mapping: Reads die meerdere matches hebben in het genoom.

  • Secondary en Duplicate: Overbodige of gedupliceerde reads, vaak door sequencing-artifacten.

  • No Features: Reads die niet aan een feature zijn gekoppeld, mogelijk door technische redenen.

We kunnen met de uitkomt van deze data op het eerste gezicht zien dat de Assigned waarden hoger zijn dan Unnasigned wat dus betekent dat de data redelijk goed is.

Visualisatie van featurecounts

Het van FeatureCounts verkregen .csv bestand kan nu als tabel in R geladen worden. De namen van de samples zijn erg lang en bevatten het gehele pad naar de bestanden, daarom worden de kolomen van de samples met de onderstaande code hernoemd naar de SSR namen.

#counts <- read.delim(file = "/home/floris/Documenten/Data_set/PRJNA885415/Transcriptomics/featurecounts/featurecounts_output.csv", header = F)

counts = read.csv( "/Users/jarnoduiker/Desktop/Transcriptomics_data/featurecounts_output.csv", row.names = "Geneid", sep="", head=T, skip=1,)

# Geef namen aan de kolomen
names(counts) <- c("Chr","Start","End","Strand","Length","SRR21754408","SRR21754417","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")

# conflict_prefer("select", "dplyr") is nodig omdat er blijkbaar twee libraries zijn die de functie "select" aanbieden.
conflict_prefer("select", "dplyr")
## [conflicted] Will prefer dplyr::select over any other package.
counts_df <- counts %>% 
  select(-c(Chr, Start, End, Strand, Length))



glimpse(counts_df)
## Rows: 41,079
## Columns: 8
## $ SRR21754408 <int> 0, 0, 1, 0, 10, 1, 0, 0, 0, 1, 438, 862, 1, 839, 10, 0, 43…
## $ SRR21754417 <int> 0, 0, 1, 1, 10, 0, 0, 0, 0, 0, 606, 1217, 3, 1244, 23, 0, …
## $ SRR21754418 <int> 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 14, 20, 0, 36, 1, 0, 9, 0, 0…
## $ SRR21754419 <int> 0, 0, 1, 0, 2, 0, 0, 0, 0, 0, 66, 127, 1, 140, 3, 0, 80, 0…
## $ SRR21754420 <int> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 44, 114, 0, 117, 0, 0, 41, 0…
## $ SRR21754421 <int> 0, 0, 1, 3, 1, 0, 1, 0, 1, 0, 61, 151, 0, 213, 5, 1, 85, 0…
## $ SRR21754422 <int> 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 96, 209, 0, 240, 1, 0, 131, …
## $ SRR21754423 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 60, 144, 0, 141, 1, 0, 94, 0…

Hier wordt de colname veranderd naar groep, en de coldata word een data frame met de aangegeven samples, (SRR21754408,SRR21754417,SRR21754418,SRR21754419,SRR21754420,SRR21754421,SRR21754422,SRR21754423) daarna worden de colnames weer aan gepast naar groep in het data frame en de kollom groep wordt ook nog factor gebruikt om een “vector” naar factor te transformen.

colnames <- c("Groep")

coldata <- data.frame(c("SRR21754408","SRR21754417","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423"), c("Jong","Jong","Jong","Jong","Oud","Oud","Oud","Oud"), row.names = 1)
colnames(coldata) <- colnames

#coldata <- coldata[,c("Groep")]

coldata$Groep <- factor(coldata$Groep)

#coldata$Type <- factor(coldata$Type)

Hier word er door het data frame geloopt en met pivot longer dit “pivot” een data frame van wijdformaat naar langformaat.

Dan wordt er ggplot gebruikt om een box plot te maken om te laten zien hoeveel reads er bij elke sample zit.

counts_long <- counts_df %>%
  pivot_longer(cols = everything(), 
               names_to = "sample_naam", 
               values_to = "reads")

ggplot(data = counts_long, aes(x = sample_naam, y = reads)) +
  geom_boxplot() + 
  xlab(" ") +
  ylab("Aantal reads:") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
Figuur 8: boxplot reads per sample

Figuur 8: boxplot reads per sample

In het plot is te zien dat het grootste deel van de samples tussen de 0 en 50000 reads zit echter is er bij samples 4419 te zien dat er 1 uitschieter is van boven de 2000000 wat wel interessant is want waarom zitten hier zoveel reads? Het kan meerdere dingen zijn van technische fouten tot biologische variant.

Echter ziet het er verder redelijk goed uit waardoor we verder kunenen

In het volgende codeblok worden de kolommen geselecteerd waarin de jonge en oude muizen staan. Deze worden in nieuwe variabellen gezet zodat we ze afhankelijk van elkaar kunnen bekijken en dus makkelijker vergelijken.

oude_muizen <- counts_df[,c(5:8)]

jonge_muizen <- counts_df[,c(1:4)]

Hier word de data van elke sample die hoort bij jong en oud in hun toebehorende data frame gezet.

# DF met alle rijen waarbij de jonge muizen een coverage hoger dan 0 hebben voor alle samples:
jonge_muizen_df <- counts[which(counts$SRR21754408 > 0 & counts$SRR21754418 > 0 & counts$SRR21754419 & counts$SRR21754417), ]

# DF met alle rijen waarbij de oude muizen een coverage hoger dan 0 hebben voor alle samples:
oude_muizen_df <- counts[which(counts$SRR21754420 > 0 & counts$SRR21754421 > 0 & counts$SRR21754422 & counts$SRR21754423 > 0), ]

In het volgende stuk code is een grafiek weergegeven van alle genen met een hogere expressie dan 0.

#ggplot word hier aangeroepen met ggplot(), dan is jonge_muizen_df als de data source gebruikt. in de aes moeten de x en y worden gespecifeerd waar de x row.names is wat dus de gen namen zijn en y de lengte is wat dus de groote is van het gen. daarna is geom_bar() aangeroepen dit geeft aan dat de plot een barplot moet zijn. xlab en ylab verzorgen de tekst op de x en y as. daarna nog het thema van hoe de grafiek er moet uitzien
ggplot(data = jonge_muizen_df, mapping = aes(x = row.names(jonge_muizen_df) , y = Length)) +
  geom_bar(stat = "identity") + 
  xlab(" ") + 
  ylab("Gene length: (bp)") +
  theme(axis.text.x = element_blank())
Figuur 9: gen expressie jonge muizen

Figuur 9: gen expressie jonge muizen

Deze plot laat de lengte van genen in baseparen (bp) zien (zoals te zien op de y as), elke balk vertegenwoordigt een gen en geeft de lengte aan. op de x-as staan de individuele genen alleen als deze weergegeven worden word de grafiek onleesbaar. Wat opvalt is dat de meeste genen relatief kort zijn vaak onder de 50.000, maar er zijn een paar genen die veel langer zijn en uitschieters vormen deze kunnen belangrijke / complexe genen zijn die vaak meer exons en introns bevatten of belangrijke regulatorische sequenties dit is dus alleen in jonge muizen.

#ggplot word hier aangeroepen met ggplot(), dan is jonge_muizen_df als de data source gebruikt. in de aes moeten de x en y worden gespecifeerd waar de x row.names is wat dus de gen namen zijn en y de lengte is wat dus de groote is van het gen. daarna is geom_bar() aangeroepen dit geeft aan dat de plot een barplot moet zijn. xlab en ylab verzorgen de tekst op de x en y as. daarna nog het thema van hoe de grafiek er moet uitzien
ggplot(data = oude_muizen_df, mapping = aes(x = row.names(oude_muizen_df) , y = Length)) +
  geom_bar(stat = "identity") + 
  xlab(" ") + 
  ylab("Gene length: (bp)") +
  theme(axis.text.x = element_blank())
Figuur 10: gen expressie oude muizen

Figuur 10: gen expressie oude muizen

Deze plot laat de lengte van genen in baseparen (bp) zien (zoals te zien op de y as), elke balk vertegenwoordigt een gen en geeft de lengte aan. op de x-as staan de individuele genen alleen als deze weergegeven worden word de grafiek onleesbaar. net zoals de vorige plot echter is wel te zien dat de lengte van de genen hier wel anders zijn dan die bij de jonge muizen dit is te zien aan de rechter kant van de x-as hier zijn 2 grote pieken te zien die niet in de jonge grafiek is te zien. Dit kan wijzen op verandering in gen lengte - sequentie lengte met leeftijd.

myColors <- hue_pal()(4)

## Plot the log2-transformed data with a 0.1 pseudocount
plotDensity(log2(counts_df + 0.1), col=rep(myColors, each=3),
            lty=c(1:ncol(counts_df)), xlab='Log2(count)',
            main='Expression Distribution')

## Add a legend and vertical line
legend('topright', names(counts_df), lty=c(1:ncol(counts_df)),
       col=rep(myColors, each=3))
abline(v=-1.5, lwd=1, col='red', lty=2)
Figuur 11: expressie distributie

Figuur 11: expressie distributie

Deze plot laat de expressie verdeling van verschillende samples (SRR21754408, SRR21754417, enzovoort) zien met log2count. Wat de log-getransformeerde read counts per gen weergeeft. op de x-as staat de genexpressie weergegeven als log2count, terwijl de y-as de dichtheid aangeeft.

de pieken rond de negatieve waarden geven waarschijlijk genen aan met lage of geen expressie in de meeste samples waar de Log2(count)-waarden positief zijn, worden genen met hogere expressie zichtbaar, hoewel dit minder frequent voorkomt.

Wat is uit te laten is dat er een grote negatieve piek is. Een negatieve log2 zegt dat er negatieve of geen expressie. De meeste genen die in onze data zitten hebben dus minder of geen expressie en zijn dus niet super interessant.

Verder is een kleine piek te zien bij 0 dit zijn dus genen die in DEZE samples niet expressed worden. als laatste is te zien dat de 4408 en 4417 sample bij de 10 ook nog expressie hebben dus bij deze samples zouden er genen kunnen zijn die overexpressed zijn.

Het aantal reads met een coverage van meer dan 0 van de jonge muizen: 14315

length(which(counts$SRR21754408 > 0 & counts$SRR21754418 > 0 & counts$SRR21754419 & counts$SRR21754417 > 0))
## [1] 14290

Het aantal reads met een coverage van meer dan 0 van de oude muizen: 15247

length(which(counts$SRR21754420 > 0 & counts$SRR21754421 > 0 & counts$SRR21754422 & counts$SRR21754423 > 0))
## [1] 15247

deseq2 Inlezen van de feature counts data inclusief de metadata collomen De deseq is gemaakt door Ivar in zijn logboek zal er uitgebreide functie uitleg zitten. DE DESEQ DATA IS NIET HELEMAAL CORRECT DIT IS OMDAT ER 1 SAMPLE MIST: SRR21754417

Deze deseq2 is gemaakt door ivar. om hier betere en verdere uitleg van te krijgen kijk naar het logboek van ivar.

# inlezen
# skip 2 lines met import data
transcrip_table <- read.delim("/Users/jarnoduiker/Desktop/Transcriptomics_data/count.out",header = FALSE,skip = 2)
# names vector van floris
names(transcrip_table) <- c("Gene_id","Chr","Start","End","Strand","Length","SRR21754408","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
#print
head(transcrip_table)
##   Gene_id
## 1 Gm26206
## 2    Xkr4
## 3 Gm53491
## 4     Rp1
## 5   Sox17
## 6 Gm22307
##                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               Chr
## 1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     NC_000067.7
## 2                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 3                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                 NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 4 NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 5                                                                                                                                                                                                                                                             NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7;NC_000067.7
## 6                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     NC_000067.7
##                                                                                                                                                                                                                                                                                                                                                                                                                                                             Start
## 1                                                                                                                                                                                                                                                                                                                                                                                                                                                         3172239
## 2                                                                                                                                                                                                                                                                                                                                                         3269956;3283662;3284705;3393039;3420932;3491925;3491925;3491925;3491925;3740775;3740775;3740775;3740775
## 3                                                                                                                                                                                                                                                                                                                                                                                                 3363896;3365064;3396936;3422903;3436480;3438011;3438565;3438565
## 4 4185896;4189891;4190089;4190238;4212835;4212835;4218035;4218035;4218835;4218835;4234078;4234078;4240428;4240428;4267757;4267757;4276883;4276883;4296834;4296834;4298666;4298666;4301276;4301276;4313356;4313356;4313640;4313640;4313766;4313766;4315254;4315254;4331750;4331750;4337692;4337692;4354989;4354989;4361069;4363149;4363149;4381493;4381493;4414373;4422133;4422133;4422133;4422133;4422425;4422425;4422425;4422425;4430423;4430423;4479393;4479393
## 5                                                                                                                                                                         4561154;4561154;4561154;4561154;4561154;4561154;4561154;4561154;4561154;4561154;4563323;4563323;4563323;4563323;4563323;4563323;4563995;4563995;4563995;4563995;4563995;4563995;4565359;4565359;4565359;4565359;4565359;4566514;4566514;4566514;4566514;4566514;4566514;4566514;4566514
## 6                                                                                                                                                                                                                                                                                                                                                                                                                                                         4599240
##                                                                                                                                                                                                                                                                                                                                                                                                                                                               End
## 1                                                                                                                                                                                                                                                                                                                                                                                                                                                         3172348
## 2                                                                                                                                                                                                                                                                                                                                                         3277540;3287191;3287191;3393983;3422536;3492124;3492124;3492124;3492124;3741733;3741733;3741733;3741721
## 3                                                                                                                                                                                                                                                                                                                                                                                                 3364050;3365091;3397056;3438060;3436728;3438060;3448035;3448035
## 4 4186171;4189935;4190296;4190296;4212989;4212989;4218186;4218186;4218967;4218967;4234164;4234164;4240627;4240627;4267864;4267864;4277060;4277060;4297046;4297046;4298842;4298842;4301367;4301367;4313485;4313485;4313671;4313671;4313842;4313842;4315329;4315329;4331828;4331828;4337843;4337843;4355121;4355121;4363235;4363235;4363235;4381656;4381656;4420314;4422304;4422304;4422304;4422304;4423060;4423060;4423060;4423060;4430526;4430526;4479508;4479489
## 5                                                                                                                                                                         4562891;4562891;4562891;4562891;4562891;4562891;4562891;4562891;4562891;4562891;4563831;4565421;4563689;4563689;4563713;4565617;4564086;4564086;4564086;4564086;4564086;4564086;4565617;4566165;4566165;4565421;4565421;4566596;4566664;4566664;4566664;4566664;4566664;4567624;4566664
## 6                                                                                                                                                                                                                                                                                                                                                                                                                                                         4599346
##                                                                                                            Strand
## 1                                                                                                               +
## 2                                                                                       -;-;-;-;-;-;-;-;-;-;-;-;-
## 3                                                                                                 +;+;+;+;+;+;+;+
## 4 -;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-
## 5                                           -;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-;-
## 6                                                                                                               +
##   Length SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421
## 1    110           0           0           0           0           0
## 2  14824           0           0           1           0           0
## 3  24933           1           1           1           1           1
## 4  12004           0           1           0           0           0
## 5   5692          10          10           0           2           0
## 6    107           1           0           0           0           0
##   SRR21754422 SRR21754423 NA
## 1           0           0  0
## 2           0           1  0
## 3           1           1  0
## 4           3           0  0
## 5           1           1  0
## 6           0           1  0

De import data veranderen naar een count matrix voor de deseq functie

Er zijn hier een paar dingen die gebeuren. Er is een nieuwe variable count_matrix, hierin word een data frame gezet die ook gelijk word aangemaakt. in dit data frame zit transcrip_table id en alle samples. De names worden ook veranderd zodat het goed leesbaar is.

count_matrix <- data.frame(transcrip_table$Gene_id, transcrip_table$SRR21754408,transcrip_table$SRR21754418,transcrip_table$SRR21754419,transcrip_table$SRR21754420,transcrip_table$SRR21754421,transcrip_table$SRR21754422,transcrip_table$SRR21754423)
names(count_matrix) <- c("Gene_ID","SRR21754408","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
head(count_matrix)
##   Gene_ID SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421
## 1 Gm26206           0           0           0           0           0
## 2    Xkr4           0           0           1           0           0
## 3 Gm53491           1           1           1           1           1
## 4     Rp1           0           1           0           0           0
## 5   Sox17          10          10           0           2           0
## 6 Gm22307           1           0           0           0           0
##   SRR21754422 SRR21754423
## 1           0           0
## 2           0           1
## 3           1           1
## 4           3           0
## 5           1           1
## 6           0           1

Hier wordt weer een dataframe aangemaakt voor de nieuwe variabele conditie_matrix. hierin staat de sample, jong of oud en de PE. Daarna word het data frame getransformeerd door de sample names niet als kol naam te zetten maar als row names en de col names aan passen naar sample, group en type. verder word de data omgezet naar “factor” van int / chr

# condition data frame
conditie_matrix <- data.frame(
                              c("SRR21754408","young","PE"),
                              c("SRR21754418","young","PE"),
                              c("SRR21754419","young","PE"),
                              c("SRR21754420","old","PE"),
                              c("SRR21754421","old","PE"),
                              c("SRR21754422","old","PE"),
                              c("SRR21754423","old","PE")
                              )
transpose <- data.frame(t(conditie_matrix))
# corigeer frame voor input functie deseq2
# condities moeten factoren zijn voor input functie
conditie_matrix <- data.frame(as.factor(transpose$X1),as.factor(transpose$X2),as.factor(transpose$X3))
                              
names(conditie_matrix) <- c("sample","group","type")
head(conditie_matrix)
##        sample group type
## 1 SRR21754408 young   PE
## 2 SRR21754418 young   PE
## 3 SRR21754419 young   PE
## 4 SRR21754420   old   PE
## 5 SRR21754421   old   PE
## 6 SRR21754422   old   PE

Hier wordt het cts data frame aangemaakt waar alle data van de samples ingezet word, hierin staat eigenlijk of de gen aanwezig is in de sample

# samples in een nieuw dataframe 
cts <- data.frame(count_matrix$SRR21754408, count_matrix$SRR21754418, count_matrix$SRR21754419, count_matrix$SRR21754420, count_matrix$SRR21754421, count_matrix$SRR21754422, count_matrix$SRR21754423)
# verander colnames en rownames
# colnames moet gelijk zijn aan rownames conditie matrix
colnames(cts) <- c("SRR21754408","SRR21754418","SRR21754419","SRR21754420","SRR21754421","SRR21754422","SRR21754423")
# rownames voor de gene names
row.names(cts) <- count_matrix$Gene_ID
head(cts)
##         SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421 SRR21754422
## Gm26206           0           0           0           0           0           0
## Xkr4              0           0           1           0           0           0
## Gm53491           1           1           1           1           1           1
## Rp1               0           1           0           0           0           3
## Sox17            10          10           0           2           0           1
## Gm22307           1           0           0           0           0           0
##         SRR21754423
## Gm26206           0
## Xkr4              1
## Gm53491           1
## Rp1               0
## Sox17             1
## Gm22307           1

Al deze voorgemaakte data frames komen nu samen om het deseq object te maken.

Deseq object aanmaak

# Cts de verwerkte fetuere counts matrix
# conditie matrix is de metadata matrix voor de vergelijking
# design is de colom (datatype=factor) van de conditie matrix.
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = conditie_matrix,
                              design = ~ group)
dds
## class: DESeqDataSet 
## dim: 41079 7 
## metadata(1): version
## assays(1): counts
## rownames(41079): Gm26206 Xkr4 ... TrnT TrnP
## rowData names(0):
## colnames(7): SRR21754408 SRR21754418 ... SRR21754422 SRR21754423
## colData names(3): sample group type

deseq analyse object

deseq <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
results(deseq)
## log2 fold change (MLE): group young vs old 
## Wald test p-value: group young vs old 
## DataFrame with 41079 rows and 6 columns
##           baseMean log2FoldChange     lfcSE      stat    pvalue      padj
##          <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
## Gm26206   0.000000             NA        NA        NA        NA        NA
## Xkr4      0.986635       0.460688   2.72789  0.168881  0.865890        NA
## Gm53491   1.726109      -1.150346   1.67134 -0.688278  0.491278        NA
## Rp1       0.533061      -1.969095   3.02457 -0.651033  0.515025        NA
## Sox17     1.314470       0.533431   1.62886  0.327486  0.743300        NA
## ...            ...            ...       ...       ...       ...       ...
## ND6     1373.89862      -0.470926  0.536920 -0.877088 0.3804389  0.923757
## TrnE       3.49194       0.338739  1.133636  0.298808 0.7650869        NA
## CYTB    4964.47156      -0.586623  0.279874 -2.096027 0.0360798  0.415000
## TrnT       2.75378      -0.902012  1.144438 -0.788170 0.4305971        NA
## TrnP      91.48406      -0.436424  0.337111 -1.294602 0.1954574  0.795742

De resultaten van de deseq geeft collommen met verschillende info

baseMean - het gemiddelde van de genormaliseerde telwaarden, gedeeld door groottefactoren, genomen over alle monsters

Log2FoldChange - De verhouding van de genormaliseerde gemiddelde gen-UMI-tellingen in elke cluster/groep ten opzichte van alle andere clusters/groepen ter vergelijking

lfcSE - De standaardfoutschatting voor de log2-voudige veranderingsschatting

stat - De waarde van de test statistiek voor het gen of het transcript

pvalue - is dus de pvalue

padj - is de adjusted pvalue voor meerdere testen voor genen of transcripten

Wat kan je dan met deze dingen. Je kan hier verschillende dingen mee zoals een kegg pathway analyse, volcano plots en andere visualisaties. Het wordt dus gebruikt om te laten zien hoe de genen worden expressed en daarmee is te zien of er een verkeerd werkend gen is.

result deseq object waar gegroupeert word waar dus de resultaten van de oude muizen samples worden vergeleken met de jonge muizen samples.

# vergelijk groep
resultaat_deseq <- results(deseq, contrast = c("group", "young","old"))
resultaat_deseq
## log2 fold change (MLE): group young vs old 
## Wald test p-value: group young vs old 
## DataFrame with 41079 rows and 6 columns
##           baseMean log2FoldChange     lfcSE      stat    pvalue      padj
##          <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
## Gm26206   0.000000             NA        NA        NA        NA        NA
## Xkr4      0.986635       0.460688   2.72789  0.168881  0.865890        NA
## Gm53491   1.726109      -1.150346   1.67134 -0.688278  0.491278        NA
## Rp1       0.533061      -1.969095   3.02457 -0.651033  0.515025        NA
## Sox17     1.314470       0.533431   1.62886  0.327486  0.743300        NA
## ...            ...            ...       ...       ...       ...       ...
## ND6     1373.89862      -0.470926  0.536920 -0.877088 0.3804389  0.923757
## TrnE       3.49194       0.338739  1.133636  0.298808 0.7650869        NA
## CYTB    4964.47156      -0.586623  0.279874 -2.096027 0.0360798  0.415000
## TrnT       2.75378      -0.902012  1.144438 -0.788170 0.4305971        NA
## TrnP      91.48406      -0.436424  0.337111 -1.294602 0.1954574  0.795742

Het is moeilijk uit te lezen op deze manier dus we zetten het in een variable zodat het kan vergeleken worden met latere visualisaties.

LFC shrinkage object

lfcShrink - gebruikt om de log fold changes (LFC) van genexpressieniveaus te verkleinen of “in te schatten”. Dit komt vaak voor in analyses met DESeq2. dit doe je omdat na de standaard DEseq2-analyse sommige genen extreem hoge of lage LFC-waarden kunnen hebben. vooral als deze lage reads hebben of grote variaties. Zulke extreme waarden kunnen ruis toevoegen aan de resultaten en om dit te vermijden past lcfShrink shrinkage toe op de LFC waarden wat dus betekent dat het deze waarden dichter naar 0 toe trekt voor genen met lage betrouwbaarheid (de genen met weinig reads of hoge variantie) dit maakt de resultaten stabieler en helpt focus op de genen met een goed-geschatte LFC-waarden

deseq_shrink <- lfcShrink(deseq, coef = "group_young_vs_old", type = "apeglm")
## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
##     Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
##     sequence count data: removing the noise and preserving large differences.
##     Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
deseq_shrink
## log2 fold change (MAP): group young vs old 
## Wald test p-value: group young vs old 
## DataFrame with 41079 rows and 5 columns
##           baseMean log2FoldChange     lfcSE    pvalue      padj
##          <numeric>      <numeric> <numeric> <numeric> <numeric>
## Gm26206   0.000000             NA        NA        NA        NA
## Xkr4      0.986635     0.00282996  0.127617  0.865890        NA
## Gm53491   1.726109    -0.00384646  0.127515  0.491278        NA
## Rp1       0.533061    -0.00356699  0.127684  0.515025        NA
## Sox17     1.314470     0.00306657  0.127293  0.743300        NA
## ...            ...            ...       ...       ...       ...
## ND6     1373.89862    -0.02464090  0.127761 0.3804389  0.923757
## TrnE       3.49194     0.00457367  0.126905 0.7650869        NA
## CYTB    4964.47156    -0.14708214  0.227104 0.0360798  0.415000
## TrnT       2.75378    -0.01102058  0.127642 0.4305971        NA
## TrnP      91.48406    -0.05977954  0.137425 0.1954574  0.795742

hoeveelheid p waardes onder de 0.01 van de 41079

hier uit is te halen dat de 43 getelde sequenties zeer hoge statistische significantie hebben.

sum(resultaat_deseq$padj < 0.01, na.rm=TRUE)
## [1] 43

Resultaat p waardes onder de 0.05

Hier worden de resultaten van deseq gefilterd op een p-waarde van < 0.05 zodat deze later kunnen worden gebruikt.

deseq_05 <- results(deseq, alpha=0.05)
summary(deseq_05)
## 
## out of 28289 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 47, 0.17%
## LFC < 0 (down)     : 138, 0.49%
## outliers [1]       : 7, 0.025%
## low counts [2]     : 16062, 57%
## (mean count < 9)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

ma plot pre shrinkage

plotMA(resultaat_deseq, ylim=c(-2,2))
Figuur 12: MA plot eerste deseq zonder shrinkage

Figuur 12: MA plot eerste deseq zonder shrinkage

ma plot met shrinkage

plotMA(deseq_shrink, ylim=c(-3,3))
Figuur 13: Ma plot eerste deseq met shrinkage

Figuur 13: Ma plot eerste deseq met shrinkage

Een MA-plot is een plot die de verhouding tussen genexpressie in twee condities laat zien. Het helpt om te identificeren welke genen significant differentieel tot expressie komen. De genen die gekleurd zijn interessant, dit is omdat ze significante genexpressie hebben dus meer expressie dan de norm (de grijze bolletjes en lijn)

volcano plot

Een volcano plot is een visualisatietool die genen laat zien die significant en sterk differentieel tot expressie komen tussen twee condities. De grootte van het expressieverschil (fold change) en de significantieniveau (padj) van elk gen, waardoor je snel kunt zien welke genen biologisch relevant zijn.

EnhancedVolcano(deseq_shrink, x="log2FoldChange", y = "padj",lab = rownames(deseq_shrink), title = "young versus old", pCutoff = 0.05,FCcutoff = 2, ylim = c(0,50),legendPosition = "none",subtitle = "FDR <= 0.05 and absolute FC >= 2")
Figuur 14:Volcano plot jong vs oud

Figuur 14:Volcano plot jong vs oud

Hier is onze volcano plot, er zijn een paar dingen te zien, het grootste deel van de genen zit tussen de standaard opgestelde voorwaarden voor singificantie (de 3 onderbroken lijnen) degene die er niet tussen zitten zijn rood gekleurd. Er zijn een paar namen te zien zoals Klra8 en Clnk wat niet de pax5 en ebf1 zijn die wij onderzoeken.

plotCounts(dds, "Pax5", intgroup = "group")
Figuur 15: pax5 expressie

Figuur 15: pax5 expressie

Hier is de genexpressie van de pax5 gen te zien per groep. De expressie van de pax5 gen is veel hoger gemiddeld in de jongere muizen wat verwacht is omdat jongere muizen nog veel ontwikkelen. er is wel 1 uitschieter bij de oudere muizen van 1600+ counts. Dit is interessant echter zou ik dit niet nu kunnen verklaren omdat er weinig verdere info is over de muizen behalve waar ze vandaan komen en dus kan het een uitschieter zijn zonder echte waarde.

plotCounts(dds, "Ebf1", intgroup = "group")
Figuur 16: expressie Ebf1 jong en oud

Figuur 16: expressie Ebf1 jong en oud

Hier is de gen epxressie te zien van de Ebf1 De expressie van de Ebf1 gen is veel hoger gemiddeld in de jongere muizen wat verwacht is omdat jongere muizen nog veel ontwikkelen. er is weer 1 uitschieter bij de oudere muizen van 6000-7000 counts wat waarschijnlijk dezelfde muis is als de vorige uitschieter en dit is dus niet uit te leggen op dit moment.

exctracting transformd value voor pca

vsd <- vst(dds, blind = FALSE)
head(assay(vsd))
##         SRR21754408 SRR21754418 SRR21754419 SRR21754420 SRR21754421 SRR21754422
## Gm26206    5.033995    5.033995    5.033995    5.033995    5.033995    5.033995
## Xkr4       5.033995    5.033995    5.641393    5.033995    5.033995    5.033995
## Gm53491    5.161996    5.136355    5.641393    5.330047    5.405110    5.308396
## Rp1        5.033995    5.136355    5.033995    5.033995    5.033995    5.507854
## Sox17      5.437583    5.357079    5.033995    5.451948    5.033995    5.308396
## Gm22307    5.161996    5.033995    5.033995    5.033995    5.033995    5.033995
##         SRR21754423
## Gm26206    5.033995
## Xkr4       5.287323
## Gm53491    5.287323
## Rp1        5.033995
## Sox17      5.287323
## Gm22307    5.287323

PCA plot waarbij je kunt zien dat de groepen verschillen en dat bij jong 1 outlyer is ssr 21754408 (gevonden met het gecomente pca plot hieronder)

# standaard pca
plotPCA(vsd, intgroup="group", ntop = 40000)
## using ntop=40000 top features by variance
Figuur 17: verschillen in groepen

Figuur 17: verschillen in groepen

# sample identificatie
#plotPCA(vsd, intgroup=c("group","sample"))

22/10/2024 KEGG pathway

Jarno’s versie van de KEGG pathway - problemen eind visualisatie.

library(clusterProfiler)
## 
## clusterProfiler v4.12.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/
## 
## Please cite:
## 
## T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan,
## X Fu, S Liu, X Bo, and G Yu. clusterProfiler 4.0: A universal
## enrichment tool for interpreting omics data. The Innovation. 2021,
## 2(3):100141
library(org.Mm.eg.db)
## Loading required package: AnnotationDbi
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:clusterProfiler':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## 
library(DOSE)
## DOSE v3.30.5  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/
## 
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609
#een lijst met alle gen entrez id's die wij nodig hebben is de 17681

search_kegg_organism('mmu', by='kegg_code')
##      kegg_code               scientific_name                   common_name
## 26        mmur            Microcebus murinus              gray mouse lemur
## 30         mmu                  Mus musculus                   house mouse
## 6386      mmuc Mycolicibacterium mucogenicum Mycolicibacterium mucogenicum
# hier is een eerste poging op het filteren van de data op de juiste pvalue deze was echter later veranderd omdat het gehele dataframe nodig was van de deseq resultaten
padj <- resultaat_deseq$padj
padj[is.na(padj)] <- 1
 
filtered_deseq_result <- rownames(resultaat_deseq)[padj < 0.05]
#head(filtered_deseq_result)

Hier worden de resultaten uit de deseq gefilterd. De manier dat we filteren is door een subset te maken. deze subset moet een P adjusted value hebben van minder dan 0.05. De reden dat we voor 0.05 hebben gekozen is omdat kort gezegt Fisher suggested is en ook dus de industry standaard.

subset_resultaat_deseq <- subset(resultaat_deseq, padj <= 0.05)

hier wordt de subset omgezet naar een dataframe dit is zodat er beter mee is te werken.

resultaat_deseq_datafr <- as.data.frame(subset_resultaat_deseq)

resultaat_deseq_datafr
##                   baseMean log2FoldChange     lfcSE      stat       pvalue
## Gm4956            9.943480     -3.3573315 0.8726760 -3.847168 1.194909e-04
## Gm41914          34.779063     -1.4053941 0.3906152 -3.597900 3.207974e-04
## Actr3          1337.546469     -0.6571379 0.1918505 -3.425261 6.142093e-04
## Rab7b            24.367399     -1.7213192 0.4883002 -3.525125 4.232827e-04
## Fasl              8.598747     -3.3404567 0.8029135 -4.160419 3.176642e-05
## Xcl1             21.771844     -2.6391352 0.6993227 -3.773845 1.607510e-04
## Fcgr2b          190.309764     -0.9868967 0.2799897 -3.524761 4.238651e-04
## Fcgr4            72.033447     -1.2747637 0.3755859 -3.394067 6.886284e-04
## Pcp4l1           94.026296      0.9150324 0.2572807  3.556553 3.757522e-04
## Slamf7          559.786247      0.7016975 0.2054318  3.415720 6.361365e-04
## Atf3             82.554561      1.6651756 0.4186536  3.977455 6.965692e-05
## Olfm1            27.374401     -1.7838464 0.4300855 -4.147656 3.358971e-05
## Gm39861          17.574010      2.9629182 0.6406062  4.625179 3.742750e-06
## Rasgrp1          30.901724     -2.0215204 0.5151597 -3.924066 8.706704e-05
## Gatm             91.795531     -1.5794845 0.4351162 -3.630029 2.833893e-04
## A630026N12Rik   242.928780      0.9573570 0.2432206  3.936168 8.279316e-05
## Dusp2           240.571632      1.5378327 0.3853270  3.990981 6.580046e-05
## Tgif2            33.922861      1.5549204 0.3881685  4.005787 6.181126e-05
## Samhd1          332.449227     -0.7677018 0.1847198 -4.156034 3.238191e-05
## Pltp             43.673993     -1.5016411 0.4239877 -3.541709 3.975437e-04
## Zfp971          114.371102     -0.9442563 0.2387260 -3.955398 7.640722e-05
## Sirpb1a          23.064096     -1.6362921 0.4606710 -3.551976 3.823502e-04
## Zbtb7b          130.096675     -1.1371122 0.2840976 -4.002541 6.266589e-05
## S100a8         5291.848808     -1.3551086 0.3066503 -4.419068 9.912751e-06
## Chil3          2856.832134     -0.9274261 0.2613861 -3.548108 3.880089e-04
## S1pr1            68.469060      0.9960447 0.2822329  3.529159 4.168829e-04
## Cnn3            242.082189      0.7292388 0.2107532  3.460155 5.398641e-04
## Gem              45.303741      1.7063905 0.4627709  3.687333 2.266165e-04
## Atp6v0d2         62.772306     -4.0316520 1.1061661 -3.644708 2.676958e-04
## Klf4            117.552904      1.2566871 0.3213334  3.910851 9.197158e-05
## Slc31a2         112.966731     -0.8872508 0.2513080 -3.530532 4.147250e-04
## Megf9           219.205452     -0.7698789 0.2237705 -3.440484 5.806743e-04
## Jun            1227.719677      0.9794741 0.2296260  4.265520 1.994374e-05
## Pde4b           470.847616      0.8981810 0.2020846  4.444579 8.806402e-06
## Zc3h12a         112.832177      0.8806101 0.2382704  3.695843 2.191585e-04
## Man1c1           72.395652     -1.2627813 0.3654650 -3.455273 5.497363e-04
## Arhgef10l        26.712355     -1.6605640 0.4875159 -3.406174 6.588010e-04
## Clnk             15.059805     -2.7826646 0.6056346 -4.594626 4.335253e-06
## Prom1            91.475427     -1.0718175 0.3158507 -3.393431 6.902303e-04
## Anxa3           163.928834     -1.0527586 0.2881733 -3.653214 2.589781e-04
## Hpse             49.523869     -1.9132288 0.4340713 -4.407638 1.045042e-05
## Btbd8            77.612103      1.2308148 0.3241262  3.797332 1.462616e-04
## Oasl2           101.904126     -0.9351533 0.2738361 -3.415011 6.377934e-04
## Ncf1            944.231018     -0.8962044 0.2517757 -3.559535 3.715115e-04
## Mgam            143.380236     -0.9896614 0.2448906 -4.041238 5.316971e-05
## Gimap3           56.929469     -1.3256432 0.3483551 -3.805436 1.415544e-04
## Gpnmb           374.813201     -6.5355986 1.2766037 -5.119520 3.063139e-07
## Skap2           153.783791     -0.7211823 0.2125233 -3.393427 6.902393e-04
## Gm44008          29.837504      1.5470014 0.4465587  3.464274 5.316653e-04
## Anxa4            62.430947     -1.5159452 0.4448542 -3.407735 6.550451e-04
## Rassf4          192.962404     -1.0461482 0.2759422 -3.791186 1.499294e-04
## C3ar1            43.136920     -1.5170326 0.4384632 -3.459886 5.404051e-04
## Clec4a3          74.859998     -1.5383463 0.4127501 -3.727065 1.937225e-04
## Klrb1a           10.476746     -2.4971028 0.7159015 -3.488054 4.865505e-04
## Cd69           1113.109624      1.4464319 0.2958760  4.888643 1.015336e-06
## Clec7a          257.033375     -0.9138534 0.2692232 -3.394408 6.877725e-04
## Klre1            15.172128     -2.2812471 0.6350108 -3.592454 3.275784e-04
## Klrk1            57.751450     -2.3554778 0.4835277 -4.871443 1.107862e-06
## Klrc3             8.230109     -3.2380495 0.9087374 -3.563240 3.663056e-04
## Klri2            29.069750     -2.6545131 0.5874070 -4.519036 6.212194e-06
## Klra8            19.849416     -3.6380844 0.7124949 -5.106120 3.288405e-07
## Klra7            19.080002     -2.8744625 0.6441256 -4.462581 8.097846e-06
## Pirb            315.668459     -1.0593741 0.2349893 -4.508181 6.538590e-06
## Ncr1             55.410115     -2.6790885 0.5715374 -4.687512 2.765462e-06
## Pglyrp1         604.584957     -1.2253921 0.3135151 -3.908559 9.284828e-05
## Cd177          1560.268810     -1.0378258 0.2933935 -3.537317 4.042143e-04
## Zfp36           901.578436      1.0917871 0.2217964  4.922475 8.545669e-07
## Rinl            125.749178     -0.9031073 0.2419821 -3.732125 1.898710e-04
## Ppp1r15a        238.156393      1.4704520 0.3159321  4.654330 3.250359e-06
## Svip             41.727998     -1.3100791 0.3823318 -3.426550 6.113001e-04
## Mfge8            92.586206     -1.7886042 0.4393962 -4.070596 4.689302e-05
## Anpep            58.089956     -3.4710640 0.9465854 -3.666932 2.454782e-04
## Gm15501          34.107488     -1.9321662 0.4947411 -3.905409 9.406617e-05
## Tenm4            36.558942     -1.8938937 0.5099552 -3.713843 2.041352e-04
## Bmal1           422.163771      0.7897508 0.2312758  3.414757 6.383891e-04
## Itpripl2        172.168920     -1.0845697 0.2924845 -3.708128 2.087975e-04
## Igsf6           232.099997     -1.1001642 0.2868806 -3.834921 1.256047e-04
## Apobr           107.168314     -0.9809726 0.2410814 -4.069052 4.720479e-05
## Itgam          1551.719408     -0.8667869 0.2177043 -3.981486 6.848561e-05
## Itgax            76.749721     -1.7231101 0.3981678 -4.327597 1.507446e-05
## Dock1            67.448952     -1.0748446 0.3172980 -3.387492 7.053467e-04
## Ifitm6          663.805623     -1.2980036 0.2970194 -4.370096 1.241918e-05
## Jund            496.680233      0.9995620 0.1789544  5.585567 2.329384e-08
## Iqcn             70.123601      1.3592055 0.3291542  4.129389 3.637291e-05
## Ddx39a          238.392529      0.8363592 0.2040125  4.099549 4.139562e-05
## Ier2            339.428271      0.8108983 0.1760569  4.605888 4.107104e-06
## Ces1d            11.919910     -2.4537084 0.7239551 -3.389310 7.006868e-04
## St3gal2          88.117861     -0.9221118 0.2668888 -3.455041 5.502090e-04
## Crispld2        101.640268     -1.0926381 0.3063453 -3.566688 3.615218e-04
## Mmp12            37.434911     -5.3030321 1.2456682 -4.257179 2.070229e-05
## Mmp8           1160.436670     -0.8628719 0.1996189 -4.322595 1.542045e-05
## Sorl1           930.504301     -0.7350288 0.2169068 -3.388685 7.022858e-04
## Il10ra          116.034134     -0.9577016 0.2686355 -3.565060 3.637731e-04
## Pkm             970.124244     -0.7488659 0.1867988 -4.008943 6.099105e-05
## Pik3cb          135.444540     -1.0868643 0.2913901 -3.729929 1.915341e-04
## Ltf            7428.912233     -0.9584815 0.2789290 -3.436291 5.897363e-04
## Eomes            10.729335     -2.3358781 0.6775369 -3.447603 5.655851e-04
## Ccr2            333.437859     -1.5840455 0.2627479 -6.028766 1.652159e-09
## Tnfaip3         491.482503      1.6742184 0.2357063  7.102986 1.220900e-12
## Sgk1            501.865759      0.8730832 0.2353681  3.709437 2.077204e-04
## Adamts14         13.489327     -2.2325827 0.6130459 -3.641787 2.707516e-04
## Itgb2          1031.779954     -0.8807052 0.2351735 -3.744917 1.804536e-04
## Prss57           41.800254     -1.9456070 0.5716883 -3.403265 6.658564e-04
## Plxnc1          245.692246     -0.8926896 0.2098783 -4.253367 2.105798e-05
## Dusp6           357.384285      1.0230410 0.2016763  5.072687 3.922377e-07
## Lyz2           4421.351744     -1.0468202 0.2422598 -4.321064 1.552789e-05
## Ppm1h           212.789840     -0.8921082 0.2335941 -3.819052 1.339654e-04
## Ddit3            62.795830      1.1257726 0.2947776  3.819058 1.339621e-04
## Mgl2             29.903115     -1.6902948 0.4566262 -3.701703 2.141568e-04
## Bcl6b            36.944070      1.7286809 0.4507191  3.835384 1.253682e-04
## Slfn8           119.432808     -0.9136515 0.2533734 -3.605949 3.110148e-04
## Slfn4           530.056281     -0.8637450 0.2210971 -3.906633 9.359121e-05
## Ccl5             29.173569     -2.2995818 0.6063877 -3.792263 1.492807e-04
## Ccl4             25.625065      1.9472130 0.4656220  4.181961 2.890060e-05
## Wfdc17           49.245118     -1.8723226 0.4727177 -3.960762 7.471086e-05
## Arl5c           281.773850      0.8116742 0.2379413  3.411237 6.466894e-04
## Ace              57.141523     -1.3204677 0.3288271 -4.015690 5.927210e-05
## Klf11            59.277023      1.0326072 0.2976904  3.468729 5.229264e-04
## Gm36723          50.184361     -1.4814379 0.3852840 -3.845055 1.205255e-04
## Gm22513         821.522087     -1.1307338 0.2838874 -3.983036 6.804052e-05
## Zfp36l1         446.088731      0.9943620 0.2118156  4.694470 2.672987e-06
## Actn1           274.870254     -0.8699495 0.2125806 -4.092329 4.270627e-05
## Fos             978.275952      1.0956212 0.2647691  4.138026 3.503072e-05
## Galc             82.891061     -1.0295785 0.2952672 -3.486938 4.885848e-04
## Ighv8-7          10.643617      2.6476451 0.7793734  3.397146 6.809265e-04
## Gm35730          28.701648      1.9056577 0.5120409  3.721690 1.978935e-04
## Aoah             99.104914     -1.1526141 0.3254953 -3.541109 3.984497e-04
## Cmah             75.928543     -1.3456313 0.3759665 -3.579126 3.447456e-04
## Irf4            295.482623      0.9087470 0.2551446  3.561695 3.684690e-04
## Tgfbi           220.152432     -0.8824823 0.2518766 -3.503629 4.589636e-04
## Arrdc3          370.110199      0.8076989 0.1957495  4.126187 3.688279e-05
## Mef2c          1586.707689      1.1568150 0.3376243  3.426338 6.117787e-04
## Vcan             35.738525     -1.7084435 0.3705706 -4.610305 4.020786e-06
## Gm36161          62.754075     -1.4636888 0.3493906 -4.189262 2.798628e-05
## Vcl             482.922822     -0.7055499 0.2031617 -3.472849 5.149645e-04
## Bmpr1a           54.588182     -2.4221931 0.4980287 -4.863562 1.152921e-06
## Lgals3          361.613544     -1.4776667 0.3179289 -4.647790 3.355095e-06
## Mmp14            46.021690     -1.8139646 0.4534762 -4.000132 6.330714e-05
## Cebpe           193.839160     -1.2958437 0.3669742 -3.531158 4.137445e-04
## Ctsg            239.796481     -1.5894319 0.4117473 -3.860212 1.132888e-04
## Adra1a            9.178068     -5.8173314 1.2057644 -4.824600 1.402843e-06
## Dock5           285.151083     -0.8542187 0.2258004 -3.783070 1.549059e-04
## Fyb1            286.793529     -0.7696619 0.2225232 -3.458793 5.426010e-04
## Zfp622          167.776550      0.8253282 0.2426455  3.401374 6.704813e-04
## Ly6c2           436.171079     -1.2025387 0.3196160 -3.762449 1.682578e-04
## Il2rb            83.646898     -2.1444046 0.4521050 -4.743156 2.104138e-06
## Bin2            229.940714     -1.1062883 0.3063897 -3.610723 3.053443e-04
## Mefv             55.631663     -1.5746766 0.3831764 -4.109534 3.964585e-05
## Ciita            72.183446     -1.4018572 0.3918927 -3.577146 3.473665e-04
## Septin5         104.380152     -1.1784543 0.3349593 -3.518201 4.344829e-04
## Nfkbiz          427.098545      1.0680589 0.3037954  3.515717 4.385677e-04
## St3gal6          97.800821      0.9761665 0.2642685  3.693844 2.208895e-04
## Fpr2            170.810100     -1.0435345 0.2982143 -3.499278 4.665207e-04
## Rpl3l            32.655640      1.5290490 0.3851948  3.969548 7.200924e-05
## Prss34          165.843413      0.8643663 0.2192175  3.942962 8.048149e-05
## Dusp1           642.988058      0.8287506 0.2445521  3.388852 7.018600e-04
## H2-Aa           450.031532     -1.0939917 0.2797218 -3.910999 9.191526e-05
## H2-Eb1          466.341001     -1.1810818 0.2868497 -4.117424 3.831312e-05
## Pla2g7          116.415223     -1.2289662 0.3253515 -3.777349 1.585063e-04
## Treml2          236.649470     -0.6263603 0.1849694 -3.386291 7.084414e-04
## Tgif1           186.413028      0.8255509 0.2167275  3.809166 1.394363e-04
## Xdh             383.203840     -0.7620958 0.2039986 -3.735790 1.871268e-04
## Egr1            365.603106      1.7818740 0.3806285  4.681399 2.849231e-06
## Cd74           1305.077804     -1.1781185 0.2951437 -3.991677 6.560766e-05
## Ahnak           811.893401     -0.9999559 0.2767688 -3.612965 3.027153e-04
## Ms4a4c           42.239897     -1.7413902 0.3711365 -4.692048 2.704832e-06
## Ms4a4b           35.218239     -2.0047348 0.5721415 -3.503914 4.584726e-04
## Ms4a3           131.068001     -1.4609811 0.4209798 -3.470431 5.196242e-04
## Mpeg1           948.120915     -1.1322084 0.3059168 -3.701034 2.147224e-04
## Anxa1          1996.122060     -0.9565209 0.2339847 -4.087964 4.351761e-05
## Gda             653.009076     -0.9597831 0.2301804 -4.169700 3.050005e-05
## Myof             84.908920     -2.1194035 0.4046852 -5.237166 1.630613e-07
## Cybb           2218.810005     -1.0088288 0.2575414 -3.917152 8.960142e-05
## G6pdx           566.541587     -1.1233134 0.2599491 -4.321283 1.551247e-05
## Mir223hg        913.079303     -1.2106112 0.3000651 -4.034495 5.471981e-05
## Vsig4            12.265848     -2.7752944 0.7360879 -3.770330 1.630318e-04
## Rnf128           13.353695     -2.9442851 0.7892219 -3.730618 1.910108e-04
## Tmsb4x         2537.406467     -0.8949135 0.2543691 -3.518169 4.345350e-04
## Tlr8             87.760130     -1.0350131 0.2914918 -3.550745 3.841425e-04
## COX1          16567.633211     -0.8470856 0.2473184 -3.425081 6.146164e-04
## COX3             16.584291     -2.0395764 0.5112723 -3.989218 6.629152e-05
##                       padj
## Gm4956        1.897904e-02
## Gm41914       3.438463e-02
## Actr3         4.809468e-02
## Rab7b         3.889496e-02
## Fasl          9.834078e-03
## Xcl1          2.204708e-02
## Fcgr2b        3.889496e-02
## Fcgr4         4.992359e-02
## Pcp4l1        3.744312e-02
## Slamf7        4.905213e-02
## Atf3          1.346173e-02
## Olfm1         9.963647e-03
## Gm39861       2.381187e-03
## Rasgrp1       1.538223e-02
## Gatm          3.143157e-02
## A630026N12Rik 1.487362e-02
## Dusp2         1.342140e-02
## Tgif2         1.342140e-02
## Samhd1        9.834078e-03
## Pltp          3.821222e-02
## Zfp971        1.412426e-02
## Sirpb1a       3.769029e-02
## Zbtb7b        1.342140e-02
## S100a8        4.515612e-03
## Chil3         3.777903e-02
## S1pr1         3.881271e-02
## Cnn3          4.494075e-02
## Gem           2.627721e-02
## Atp6v0d2      3.021646e-02
## Klf4          1.538223e-02
## Slc31a2       3.881271e-02
## Megf9         4.687659e-02
## Jun           7.259311e-03
## Pde4b         4.160209e-03
## Zc3h12a       2.584813e-02
## Man1c1        4.498664e-02
## Arhgef10l     4.972193e-02
## Clnk          2.404181e-03
## Prom1         4.992359e-02
## Anxa3         2.949344e-02
## Hpse          4.596381e-03
## Btbd8         2.119962e-02
## Oasl2         4.905213e-02
## Ncf1          3.731204e-02
## Mgam          1.255888e-02
## Gimap3        2.075318e-02
## Gpnmb         6.990600e-04
## Skap2         4.992359e-02
## Gm44008       4.490987e-02
## Anxa4         4.972193e-02
## Rassf4        2.124832e-02
## C3ar1         4.494075e-02
## Clec4a3       2.446466e-02
## Klrb1a        4.239387e-02
## Cd69          1.336864e-03
## Clec7a        4.992359e-02
## Klre1         3.481885e-02
## Klrk1         1.336864e-03
## Klrc3         3.730017e-02
## Klri2         3.301522e-03
## Klra8         6.990600e-04
## Klra7         3.972616e-03
## Pirb          3.335989e-03
## Ncr1          2.137761e-03
## Pglyrp1       1.538223e-02
## Cd177         3.847578e-02
## Zfp36         1.336864e-03
## Rinl          2.443017e-02
## Ppp1r15a      2.252328e-03
## Svip          4.809468e-02
## Mfge8         1.136032e-02
## Anpep         2.820788e-02
## Gm15501       1.538223e-02
## Tenm4         2.527908e-02
## Bmal1         4.905213e-02
## Itpripl2      2.536392e-02
## Igsf6         1.930227e-02
## Apobr         1.136032e-02
## Itgam         1.343898e-02
## Itgax         5.825243e-03
## Dock1         4.992359e-02
## Ifitm6        5.280221e-03
## Jund          9.903762e-05
## Iqcn          1.022696e-02
## Ddx39a        1.077553e-02
## Ier2          2.381187e-03
## Ces1d         4.992359e-02
## St3gal2       4.498664e-02
## Crispld2      3.730017e-02
## Mmp12         7.259311e-03
## Mmp8          5.825243e-03
## Sorl1         4.992359e-02
## Il10ra        3.730017e-02
## Pkm           1.342140e-02
## Pik3cb        2.443017e-02
## Ltf           4.730872e-02
## Eomes         4.594929e-02
## Ccr2          1.053664e-05
## Tnfaip3       1.557257e-08
## Sgk1          2.536392e-02
## Adamts14      3.029330e-02
## Itgb2         2.397589e-02
## Prss57        4.992359e-02
## Plxnc1        7.259311e-03
## Dusp6         7.147130e-04
## Lyz2          5.825243e-03
## Ppm1h         2.010270e-02
## Ddit3         2.010270e-02
## Mgl2          2.559611e-02
## Bcl6b         1.930227e-02
## Slfn8         3.361859e-02
## Slfn4         1.538223e-02
## Ccl5          2.124832e-02
## Ccl4          9.451977e-03
## Wfdc17        1.401378e-02
## Arl5c         4.939235e-02
## Ace           1.342140e-02
## Klf11         4.446617e-02
## Gm36723       1.897904e-02
## Gm22513       1.343898e-02
## Zfp36l1       2.137761e-03
## Actn1         1.088367e-02
## Fos           1.015493e-02
## Galc          4.239387e-02
## Ighv8-7       4.992359e-02
## Gm35730       2.474639e-02
## Aoah          3.821222e-02
## Cmah          3.631688e-02
## Irf4          3.730017e-02
## Tgfbi         4.065334e-02
## Arrdc3        1.022696e-02
## Mef2c         4.809468e-02
## Vcan          2.381187e-03
## Gm36161       9.393817e-03
## Vcl           4.438090e-02
## Bmpr1a        1.336864e-03
## Lgals3        2.252328e-03
## Mmp14         1.342140e-02
## Cebpe         3.881271e-02
## Ctsg          1.829111e-02
## Adra1a        1.491106e-03
## Dock5         2.171235e-02
## Fyb1          4.494075e-02
## Zfp622        4.992359e-02
## Ly6c2         2.259083e-02
## Il2rb         2.064484e-03
## Bin2          3.328774e-02
## Mefv          1.053506e-02
## Ciita         3.631688e-02
## Septin5       3.930847e-02
## Nfkbiz        3.939388e-02
## St3gal6       2.584813e-02
## Fpr2          4.103774e-02
## Rpl3l         1.370863e-02
## Prss34        1.466488e-02
## Dusp1         4.992359e-02
## H2-Aa         1.538223e-02
## H2-Eb1        1.039753e-02
## Pla2g7        2.197552e-02
## Treml2        4.992359e-02
## Tgif1         2.068035e-02
## Xdh           2.443017e-02
## Egr1          2.137761e-03
## Cd74          1.342140e-02
## Ahnak         3.328563e-02
## Ms4a4c        2.137761e-03
## Ms4a4b        4.065334e-02
## Ms4a3         4.446617e-02
## Mpeg1         2.559611e-02
## Anxa1         1.088367e-02
## Gda           9.725703e-03
## Myof          5.199617e-04
## Cybb          1.538223e-02
## G6pdx         5.825243e-03
## Mir223hg      1.269002e-02
## Vsig4         2.212202e-02
## Rnf128        2.443017e-02
## Tmsb4x        3.930847e-02
## Tlr8          3.769029e-02
## COX1          4.809468e-02
## COX3          1.342140e-02

Hier beginnen we met gebruiken van KEGG. KEGG ook wel Kyoto Encyclopedia of Genes and Genomes is een website die door middel van clusterprofiler gebruikt kan worden in Rstudio. clusterprofiler gebruikt de KEGG database voor info en “connect” met hun server om informatie te verkrijgen. Dit is te zien door een analyse te runnen en in de terminal is te zien dat er communicatie is tussen de server en jouw R.

Hier onder zoek ik de gene symbols uit. De gene symbols zijn nodig om entrez_id’s te vinden. De entrez id’s zijn nodig in latere stappen en deze entrez_id’s zijn te linken aan de gene symbolen.

#hier maak ik een data frame aan waar ik de gen namen in zet en de pvalues ook overneem
data_kegg <- data.frame(
  #De gen naam is nodig voor het verkrijgen van de ENTREZID
  Gene = row.names(deseq_shrink),
  #Pvalue is nodig voor filtering in volgende stappen 
  pvalue = deseq_shrink$pvalue,
  stringsAsFactors = FALSE  
)


gene_symbols <- row.names(deseq_shrink)

head(gene_symbols)
## [1] "Gm26206" "Xkr4"    "Gm53491" "Rp1"     "Sox17"   "Gm22307"

Met bitr kunnen we doormiddel van fromType en toType de gene symbolen omzetten naar de gewilde entrez_id’s verder is nog de Org.Mm.eg.db dit betekent Organism Mouse musculus eg database. We willen natuurlijk de entrez id’s van de muis terug krijgen en niet de mens daarom specificeren we de organism database. De entrez id’s zijn essentieel voor de KEGG analyse.

#hier gebruik ik bitr dit is een Biological Id TRanslator, hiermee vertaal ik mijn gen namen naar de ENTREID die de kegg tools verwachten met de data base van de mmu
entrez_ids <- bitr(rownames(resultaat_deseq_datafr), fromType = "SYMBOL", 
                   toType = "ENTREZID", 
                   OrgDb = org.Mm.eg.db)
## 'select()' returned 1:1 mapping between keys and columns
head(entrez_ids)
##    SYMBOL  ENTREZID
## 1  Gm4956    241041
## 2 Gm41914 105246657
## 3   Actr3     74117
## 4   Rab7b    226421
## 5    Fasl     14103
## 6    Xcl1     16963
#entrez_ids$ENTREID geeft aan de gene de juiste id's mee voor elk gen, orgbd vraagt welke database org.Mm.eg.db is de database van de mus musculus, keytipe geeft aan dat de genen die worden gebruikt in de analyse geïdentificeerd met entrez ID's. ont = BP geeft aan dat wel Biologische processen willen. Door dit te veranderen kan je ook MF - moleculaire functie krijgen en CC - cellular component. pAdjustMethod geeft aan dat de Benjamini-Hochberg (BH) methode wordt gebruikt om de p-waarden voor meerdere tests te corrigeren. Deze methode helpt om het risico op vals-positieve resultaten te verminderen. qvalueCutoff = 0.05 stelt dat alles onder 0.05 statistisch significant is en dus belangrijk is. readable = true zet de gen symbolen inplaats van de entrez IDs zodat het te bekijken is.
go_results <- enrichGO(gene = entrez_ids$ENTREZID,
                        OrgDb = org.Mm.eg.db,
                        keyType = "ENTREZID",
                        ont = "BP",
                        pAdjustMethod = "BH",
                        qvalueCutoff = 0.05,
                        readable = TRUE) 

head(go_results)
##                    ID                           Description GeneRatio   BgRatio
## GO:0002697 GO:0002697 regulation of immune effector process    25/168 473/28905
## GO:0045088 GO:0045088  regulation of innate immune response    24/168 486/28905
## GO:0097529 GO:0097529           myeloid leukocyte migration    18/168 255/28905
## GO:0007159 GO:0007159          leukocyte cell-cell adhesion    22/168 445/28905
## GO:0002274 GO:0002274          myeloid leukocyte activation    18/168 281/28905
## GO:0030595 GO:0030595                  leukocyte chemotaxis    17/168 242/28905
##            RichFactor FoldEnrichment   zScore       pvalue     p.adjust
## GO:0002697 0.05285412       9.093741 13.57029 5.797906e-17 1.784595e-13
## GO:0045088 0.04938272       8.496473 12.74335 1.165385e-15 1.793528e-12
## GO:0097529 0.07058824      12.144958 13.66782 1.219735e-14 1.251448e-11
## GO:0007159 0.04943820       8.506019 12.20071 1.921520e-14 1.478610e-11
## GO:0002274 0.06405694      11.021225 12.90689 6.524299e-14 3.816602e-11
## GO:0030595 0.07024793      12.086408 13.24191 7.439769e-14 3.816602e-11
##                  qvalue
## GO:0002697 1.171177e-13
## GO:0045088 1.177039e-12
## GO:0097529 8.212885e-12
## GO:0007159 9.703677e-12
## GO:0002274 2.504722e-11
## GO:0030595 2.504722e-11
##                                                                                                                                                                     geneID
## GO:0002697 Xcl1/Fcgr2b/Rasgrp1/Zbtb7b/Zc3h12a/Clnk/Ncf1/Gimap3/Klrb1a/Clec7a/Klre1/Klrk1/Klrc3/Klri2/Pglyrp1/Cd177/Itgam/Ccr2/Tnfaip3/Itgb2/Lgals3/Nfkbiz/Cd74/Anxa1/Vsig4
## GO:0045088                 Rab7b/Xcl1/Rasgrp1/Samhd1/S100a8/Clnk/Ncf1/Gimap3/Klrb1a/Clec7a/Klre1/Klrk1/Klrc3/Klri2/Mmp12/Ltf/Tnfaip3/Irf4/Mefv/Nfkbiz/Fpr2/Cd74/Vsig4/Tlr8
## GO:0097529                                                         Xcl1/S100a8/Pde4b/C3ar1/Cd177/Itgam/Ccr2/Itgb2/Ccl5/Ccl4/Lgals3/Mmp14/Ctsg/Fpr2/Dusp1/Pla2g7/Cd74/Anxa1
## GO:0007159                       Xcl1/Rasgrp1/Sirpb1a/Zbtb7b/S100a8/Klf4/Zc3h12a/Gimap3/Gpnmb/Cd177/Itgam/Ccr2/Itgb2/Ccl5/Lgals3/Ctsg/Nfkbiz/H2-Aa/H2-Eb1/Cd74/Anxa1/Vsig4
## GO:0002274                                                            Fcgr4/Rasgrp1/Jun/Clnk/Anxa3/Gimap3/Klrk1/Pirb/Cd177/Itgam/Jund/Mmp8/Ccr2/Itgb2/Ccl5/Irf4/Ctsg/Anxa1
## GO:0030595                                                               Xcl1/S100a8/S1pr1/Pde4b/C3ar1/Itgam/Ccr2/Itgb2/Ccl5/Ccl4/Lgals3/Ctsg/Fpr2/Dusp1/Pla2g7/Cd74/Anxa1
##            Count
## GO:0002697    25
## GO:0045088    24
## GO:0097529    18
## GO:0007159    22
## GO:0002274    18
## GO:0030595    17

De GO-verrijkingsanalyse helpt te begrijpen welke biologische processen, functies of cellulaire componenten oververtegenwoordigd zijn in een bepaalde genenlijst. Dit geeft inzicht in biologische context van de genen die bestuurt zijn

Het uitvoeren van een GO-verrijkingsanalyse helpt onderzoekers te begrijpen welke biologische processen, functies of cellulaire componenten oververtegenwoordigd zijn in een bepaalde genenlijst. Dit kan waardevolle inzichten opleveren over de biologische context van de genen die je hebt bestudeerd.

De volgende data wordt terug gegeven GeneRatio - het percentage van het aantal genen dat aanwezig is in deze GO-term ten opzichte van het totale aantal genen in deze categorie.

BgRatio - verhouding van alle genen die in deze term zijn geannoteerd

RichFactor - De verhouding vanhet differentieel tot expressie gebrachte eiwitnummer dat in deze padterm is geannoteerd tot alle eiwitnummers die zijn geannoteerd

FoldEnrichment - het percentage genen in uw lijst dat tot een pad behoort, gedeeld door het overeenkomstige percentage in de achtergrond

zScore - geeft het aantal standaarddeviaties weer dat een gegevenspunt van het gemiddelde verwijderd is als een standaardscore ID - geeft het Gene Ontology ID weer wat kan worden opgezocht om te zien wat dit gen doet - bijvoorbeeld GO:0051051 in google geeft een website met uitleg wat het doet voor functie : negative regulation of transport https://www.informatics.jax.org/vocab/gene_ontology/GO:0051051

Description - geeft ook weer wat er te vinden is over de gen functie

geneID - geeft welke genen bij deze functie werken.

#entrez_ids$ENTREID geeft aan de gene de juiste id's mee voor elk gen, orgbd vraagt welke database org.Mm.eg.db is de database van de mus musculus,
ego <- enrichGO(gene          = entrez_ids$ENTREZID,
                universe      = names(entrez_ids$SYMBOL),
                OrgDb         = org.Mm.eg.db,
                ont           = "CC",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.01,
                qvalueCutoff  = 0.05,
        readable      = TRUE)
head(ego)
##                    ID                  Description GeneRatio   BgRatio
## GO:0002102 GO:0002102                     podosome     4/170  30/28618
## GO:0070062 GO:0070062        extracellular exosome     6/170 106/28618
## GO:0042613 GO:0042613 MHC class II protein complex     3/170  13/28618
## GO:0045121 GO:0045121                membrane raft     9/170 305/28618
## GO:0098857 GO:0098857         membrane microdomain     9/170 306/28618
## GO:1903561 GO:1903561        extracellular vesicle     6/170 128/28618
##            RichFactor FoldEnrichment    zScore       pvalue    p.adjust
## GO:0002102 0.13333333      22.445490  9.084800 2.919493e-05 0.004417486
## GO:0070062 0.05660377       9.528746  6.800407 4.199625e-05 0.004417486
## GO:0042613 0.23076923      38.847964 10.551249 5.637809e-05 0.004417486
## GO:0045121 0.02950820       4.967445  5.384910 9.291897e-05 0.004417486
## GO:0098857 0.02941176       4.951211  5.371756 9.525372e-05 0.004417486
## GO:1903561 0.04687500       7.890993  6.040198 1.199213e-04 0.004417486
##                 qvalue                                            geneID Count
## GO:0002102 0.003368494                              Actr3/Vcl/Dock5/Bin2     4
## GO:0070062 0.003368494                  Fasl/Prom1/Anpep/Ace/Ahnak/Anxa1     6
## GO:0042613 0.003368494                                 H2-Aa/H2-Eb1/Cd74     3
## GO:0045121 0.003368494 Fasl/Hpse/Cd177/Itgam/Itgb2/Vcl/Bmpr1a/Ahnak/Myof     9
## GO:0098857 0.003368494 Fasl/Hpse/Cd177/Itgam/Itgb2/Vcl/Bmpr1a/Ahnak/Myof     9
## GO:1903561 0.003368494                  Fasl/Prom1/Anpep/Ace/Ahnak/Anxa1     6

Hier word de enrichKEGG gemaakt. De enrich kegg zoekt uit welke genen statistich significant zijn in hun expressie, er zijn 3 dingen nodig om de analyse te doen. Er moeten entrez_ids zijn deze zijn gemaakt in een eerdere functie, organisme moet worden gespecificeerd en de pvalueCutoff moet worden aangegeven die weer 0.05 is.

De resultaten zet ik in de variabele kk en deze zet ik om in een dataframe zodat deze makkelijker te lezen is.

kk <- enrichKEGG(gene         = entrez_ids$ENTREZID,
                 organism     = 'mmu',
                 pvalueCutoff = 0.05)
## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/mmu"...
kk
## #
## # over-representation test
## #
## #...@organism     mmu 
## #...@ontology     KEGG 
## #...@keytype      kegg 
## #...@gene     chr [1:181] "241041" "105246657" "74117" "226421" "14103" "16963" "14130" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...38 enriched terms found
## 'data.frame':    38 obs. of  14 variables:
##  $ category      : chr  "Human Diseases" "Human Diseases" "Cellular Processes" "Organismal Systems" ...
##  $ subcategory   : chr  "Infectious disease: parasitic" "Infectious disease: bacterial" "Transport and catabolism" "Development and regeneration" ...
##  $ ID            : chr  "mmu05140" "mmu05152" "mmu04145" "mmu04380" ...
##  $ Description   : chr  "Leishmaniasis - Mus musculus (house mouse)" "Tuberculosis - Mus musculus (house mouse)" "Phagosome - Mus musculus (house mouse)" "Osteoclast differentiation - Mus musculus (house mouse)" ...
##  $ GeneRatio     : chr  "9/103" "12/103" "11/103" "9/103" ...
##  $ BgRatio       : chr  "70/9793" "180/9793" "183/9793" "136/9793" ...
##  $ RichFactor    : num  0.1286 0.0667 0.0601 0.0662 0.0805 ...
##  $ FoldEnrichment: num  12.22 6.34 5.72 6.29 7.65 ...
##  $ zScore        : num  9.72 7.45 6.64 6.41 6.42 ...
##  $ pvalue        : num  4.22e-08 3.68e-07 3.22e-06 1.23e-05 3.40e-05 ...
##  $ p.adjust      : num  9.04e-06 3.93e-05 2.30e-04 6.56e-04 1.46e-03 ...
##  $ qvalue        : num  6.80e-06 2.96e-05 1.73e-04 4.94e-04 1.10e-03 ...
##  $ geneID        : chr  "246256/16476/17969/16409/16414/14281/14960/14969/13058" "14130/246256/242341/56644/16409/16411/16154/16414/12265/14960/14969/16149" "226421/14130/246256/242341/17969/56644/16409/16414/14960/14969/13058" "14130/246256/320832/16476/17969/18733/16478/74769/14281" ...
##  $ Count         : int  9 12 11 9 7 8 8 10 8 7 ...
## #...Citation
##  T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
##  clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
##  The Innovation. 2021, 2(3):100141
kk_data_frame <- as.data.frame(kk)

hier is de module over-representation analysis kegg analyse te vinden. Met de pvalueCutoff en de qvalueCutoff

Er zijn verschillende dingen te zien in deze dataframe

GeneRatio - het percentage van het aantal genen dat aanwezig is in deze GO-term ten opzichte van het totale aantal genen in deze categorie.

Bgratio - verhouding van alle genen die in deze term zijn geannoteerd

RichFactor - De verhouding vanhet differentieel tot expressie gebrachte eiwitnummer dat in deze padterm is geannoteerd tot alle eiwitnummers die zijn geannoteerd

FoldEnrichment - het percentage genen in uw lijst dat tot een pad behoort, gedeeld door het overeenkomstige percentage in de achtergrond

zScore - geeft aan hoeveel standaarddeviaties een score van het gemiddelde af zit

pvalue - is een getal tussen 0 en 1, waarmee je bepaalt of een steekproefuitkomst statistisch significant is.

geneID - geeft welke genen bij deze functie werken.

count - hoeveel er aanwezig zijn die voldoen aan de paramaters

mkk <- enrichMKEGG(gene = entrez_ids$ENTREZID,
                   organism = 'mmu',
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.05)
## Reading KEGG annotation online: "https://rest.kegg.jp/link/mmu/module"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/module"...
mkk_data_frame <- as.data.frame(mkk)

mkk_data_frame
##            ID                                      Description GeneRatio
## M00959 M00959 Guanine ribonucleotide degradation, GMP => Urate      2/10
##        BgRatio RichFactor FoldEnrichment   zScore      pvalue   p.adjust
## M00959  12/886  0.1666667       14.76667 5.127236 0.007130478 0.04991335
##            qvalue      geneID Count
## M00959 0.03752883 22436/14544     2

Er is te zien dat er maar 1 gen uit de enrichMKEGG komt en dit is de M00959, deze gen heeft te maken met de Huanine ribonucleotide degredation. Deze gen is niet vernoemd in ons onderzoek waar wij ons onderzoek op baseren en om deze reden zullen wij niks verder doen met de data die hieruit is gekomen.

23-10-2024

#hier zet ik een lege lijst in de kk_list_id
kk_list_id <- list()

#hier neem ik alle ID's die uit de kegg analyse komen die belangrijk zijn en zet deze in de nieuwe lijst
kk_list_id <- kk_data_frame$ID

kk_list_id
##  [1] "mmu05140" "mmu05152" "mmu04145" "mmu04380" "mmu05323" "mmu04650"
##  [7] "mmu05150" "mmu04613" "mmu05135" "mmu05146" "mmu04928" "mmu04670"
## [13] "mmu04666" "mmu04061" "mmu04810" "mmu05142" "mmu04620" "mmu04659"
## [19] "mmu05162" "mmu05418" "mmu05166" "mmu04010" "mmu04668" "mmu04932"
## [25] "mmu04148" "mmu04662" "mmu05417" "mmu04658" "mmu05202" "mmu04657"
## [31] "mmu04610" "mmu05221" "mmu04614" "mmu05415" "mmu05020" "mmu05171"
## [37] "mmu05100" "mmu05133"

De reden om deze lijst te maken is het makkelijk bekijken van welke gen pathways interresant zijn.

De uiteindlijk interresante in mijn analyse zijn de volgende pathways.

#mmu04666: B-celreceptoren
#mmu04612: B-celactivatie
#mmu04611: B-celdifferentiatie
#mmu04670: B-celontwikkeling

#hier word de bijbehorende pathway geopend in de KEGG website.
browseKEGG(kk, "mmu04611")

Hier worden de pathway’s met de expressie data aangeroepen wat hier uitkomt zijn foto’s. Deze foto zoals genoemd in de inleiding hebben een bepaalde structuur.

het belangerijkste om te noemen is de expressie schaal wat dus de kleuren schaal is. De kleuren schaal en cijfers betekenen ook wat, wanneer de gen een witte/grijze kleur heeft is de verandering in de expressie 0 wat dus betekent dat de expressie van dit gen normaal is. Wanneer een gen block een groene kleur heeft en dus richting de -1 gaat zegt dit dat de expressie is minder is dan de norm. Wanneer een gen een rode kleur heeft en dus richting de 1 gaat zegt dit dat het een overexpressie heeft meer dan het normaal.

#hierw worden alle pathviews 
library(pathview)
## 
## ##############################################################################
## Pathview is an open source software package distributed under GNU General
## Public License version 3 (GPLv3). Details of GPLv3 is available at
## http://www.gnu.org/licenses/gpl-3.0.html. Particullary, users are required to
## formally cite the original Pathview paper (not just mention it) in publications
## or products. For details, do citation("pathview") within R.
## 
## The pathview downloads and uses KEGG data. Non-academic uses may require a KEGG
## license agreement (details at http://www.kegg.jp/kegg/legal.html).
## ##############################################################################
mmu04612_pathview <- pathview(gene.data  = resultaat_deseq_datafr,
                     pathway.id = "mmu04612",
                     species    = "mmu")
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04612.pathview.png
mmu04666_pathview <- pathview(gene.data  = resultaat_deseq_datafr,
                     pathway.id = "mmu04666",
                     species    = "mmu")
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04666.pathview.png
mmu04611_pathview <- pathview(gene.data  = resultaat_deseq_datafr,
                     pathway.id = "mmu04611",
                     species    = "mmu")
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04611.pathview.png
mmu04670_pathview <- pathview(gene.data  = resultaat_deseq_datafr,
                     pathway.id = "mmu04670",
                     species    = "mmu")
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04670.pathview.png

De visualisaties zijn te vinden in de github map, deze hebben dus incorrecte schaal en de expressie hierin is dus niet goed.

Floris zijn versie van de KEGG pathway - problemen met visualisatie opgelost.

Deseq - floris

Statistische bewerkingen die door DESeq2 uitgevoerd worden:

Een van de tests zorgt er voor dat er rekening gehouden wordt met de false discovery rate (FDR). Hiervoor is de Benjamini Hockberg compensatie uitgevoerd om hier rekening mee te houden. De false discovery rate is het aantal genen in dit geval die ondanks dat ze niet echt differentialy expressed zijn toch zo lijken. Dit zijn dus false positieven. Wanneer de p-waarde bijvoorbeeld 0.01 is en er 4000 genen zijn, zijn er 40 false positieven bij. De kans is namelijk 1% met een p-waarde van 0.01 of lager.

levels(coldata$Groep) geeft “Jong” als eerste factor, dat betekend dat up- en downregulatie tenopzichte van de jonge groep bepaald zal worden.

dds <- DESeqDataSetFromMatrix(countData = counts_df,
                              colData = coldata,
                              design= ~ Groep)

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
resultsNames(dds) # De naam is: Groep_Oud_vs_Jong
## [1] "Intercept"         "Groep_Oud_vs_Jong"
res <- results(dds, name="Groep_Oud_vs_Jong")

# Gebruikt lfcschrink om de data met Log2 te transformeren:
res <- lfcShrink(dds, coef="Groep_Oud_vs_Jong", type="apeglm")
## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
##     Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
##     sequence count data: removing the noise and preserving large differences.
##     Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
summary(res)
## 
## out of 28500 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 1587, 5.6%
## LFC < 0 (down)     : 1142, 4%
## outliers [1]       : 3, 0.011%
## low counts [2]     : 12953, 45%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

De betekenis van de kolommen van het resultaat van DESeq2:

mcols(res)$description
## [1] "mean of normalized counts for all samples"
## [2] "log2 fold change (MAP): Groep Oud vs Jong"
## [3] "posterior SD: Groep Oud vs Jong"          
## [4] "Wald test p-value: Groep Oud vs Jong"     
## [5] "BH adjusted p-values"

Visualisatie van Differentially expressed genes:

MA-Plot

plotMA(res, ylim=c(-2,2))
Figuur 18: Ma plot 2de deseq

Figuur 18: Ma plot 2de deseq

Hier is een MA-plot te zien waar dus nog heel veel ruis in zit waardoor er geen goede conclusies uit zijn te halen.

plotCounts(dds, gene=which.min(res$padj), intgroup="Groep")
Figuur 19: gen met kleinste p waarde

Figuur 19: gen met kleinste p waarde

Visualisatie van het gen met de kleinste p-waarde: dat is: Gpnmb. - er is hier te zien dat oude muizen dit gen dus hoog expressen en de jonge muizen amper.

hieronder is het gen te zien waarop deze grafiek is gebaseerd.

counts_df[which.min(res$padj),]
##       SRR21754408 SRR21754417 SRR21754418 SRR21754419 SRR21754420 SRR21754421
## Gpnmb          25          50           1           6         277         816
##       SRR21754422 SRR21754423
## Gpnmb        1008         697

Plot van relevante genen:

Efb1, pax5 en foxo5

d <- plotCounts(dds, gene="Ebf1", intgroup="Groep", 
                returnData=TRUE)
ggplot(d, aes(x=Groep, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))
Figuur 20: Ebf1 counts nieuwe deseq

Figuur 20: Ebf1 counts nieuwe deseq

Hier is weer een genexpressie plot te zien van gen Ebf1. er is weer te zien dat de jonge muizen hoog in de expressie en de oude laag maar de uitschieter is weg bij de oude muizen.

d <- plotCounts(dds, gene="Pax5", intgroup="Groep", 
                returnData=TRUE)
ggplot(d, aes(x=Groep, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))
Figuur 21: Pax5 count nieuwe deseq

Figuur 21: Pax5 count nieuwe deseq

Hier is weer een genexpressie plot te zien van gen Pax5. er is weer te zien dat de jonge muizen hoog in de expressie en de oude laag maar de uitschieter is weg bij de oude muizen. verder is er een oude muis die heel lage expressie in het algemeen heeft dus deze vergelijken met de andere samples zou ook nog interresant kunnen zijn.

# Plot gevonden op de website: https://introtogenomics.readthedocs.io/en/latest/2021.11.11.DeseqTutorial.html

vsd <- vst(dds)

#head(assay(vsd), 30)

select <- order(rowMeans(counts(dds,normalized=TRUE)),
                decreasing=TRUE)[1:10]
df <- as.data.frame(colData(dds))
pheatmap(assay(vsd)[select,], annotation_col=df)
Figuur 22: heatmap genen

Figuur 22: heatmap genen

Dit is een pheatmap plot. Het is een visualisatie van gegevens in de vorm van een hittekaart, waarbij de waarden van een matrix visueel worden weergegeven met kleuren. Dit is handig om patronen, correlaties en trends in data te identificeren, vooral wanneer je te maken hebt met veel variabelen of een grote dataset.

Uit deze heatmap is te halen dat jong gesplit is bij de gm26917 gen dit is te zien door de donkere kleur bij de ene helft jonge muozen en de andere jonge muizen zijn geel. verder zou hier nog uit te halen zijn dat er 2 jonge muizen zijn met hoge expressie voor deze genen en de andere 2 jonge muizen lage gen epxressie hebben voor deze genen.

sampledists <- dist(t(assay(vsd)))

# Convert the 'dist' object into a matrix for creating a heatmap
sampleDistMatrix <- as.matrix(sampledists)
                         
pheatmap(sampleDistMatrix, show_colnames = FALSE,
         annotation_col = coldata, # Gebruikt dezelfde coldata die ook voor DESeq2 gebruikt is.
         clustering_distance_rows = sampledists,
         clustering_distance_cols = sampledists,
         main = "Euclidean Sample Distances")
Figuur 23: Euclidean sample distance

Figuur 23: Euclidean sample distance

Euclidean Sample Distances is een manier om de afstand te berekenen tussen twee punten in een ruimte. In de context van sample distances (afstanden tussen monsters of datapunten) verwijst het naar de afstand tussen twee gegevenspunten.

er is te zien dat de sample 4418 en 4420 een donker rode vlak hebben wat kan beteken dat ze ver uit elkaar liggen.

Principal component analysis:

Hier is een pcaplot te zien, een PCA plot is een grafische weergave van de resultaten van Principal Component Analysis (PCA), een techniek die wordt gebruikt om de meest significante patronen in een dataset te identificeren. De plot toont meestal de eerste twee of drie principale componenten als assen, waardoor we de data kunnen visualiseren in lagere dimensies. Dit helpt om structuren, clusters of trends in de data te ontdekken, zelfs wanneer de oorspronkelijke dataset veel variabelen bevat.

plotPCA(vsd, intgroup=c("Groep"))
## using ntop=500 top features by variance
Figuur 24: PCA

Figuur 24: PCA

In onze pcaplot is te zien dat er enorm grote variatie zit tussen de jonge samples, wel 15-20 punten.

degPlot(dds = dds, res = res, n = 6, xs = "Groep")
## No genes were mapped to rowData. check ann parameter values.
## Using gene as id variables
## `geom_smooth()` using formula = 'y ~ x'
Figuur 25: top expressed genen

Figuur 25: top expressed genen

De degPlot() function weergeeft volgens de documentatie enkele “top genen” en de expressie verschillen tussen jong en oud. opmerkelijk is dat alle genen in de oude muizen hogere expressie hebben.

deseq.volcano <- function(res, datasetName) {
  return(EnhancedVolcano(res, x = 'log2FoldChange', y = 'padj',
                         lab=rownames(res),
                         title = paste(datasetName, "Jonge vs Oude muizen"),
                         subtitle = bquote(italic('FDR <= 0.05 and absolute FC >= 2')),
                         # Change text and icon sizes
                         labSize = 3, pointSize = 1.5, axisLabSize=10, titleLabSize=12,
                         subtitleLabSize=8, captionLabSize=10,
                         # Disable legend
                         legendPosition = "none",
                         # Set cutoffs
                         pCutoff = 0.05, FCcutoff = 2))
}

deseq.volcano(res = res, datasetName = "Volcano plot:")
Figuur 26: Volcano plot: Jonge vs Oude muizen

Figuur 26: Volcano plot: Jonge vs Oude muizen

De analyse van floris

Een volcano plot weergeeft alle genen maar enkel de genen die statistisch significant zijn en een hoge fold change hebben worden rood gekleurd. De genen die wel statistisch significant zijn maar geen significante fold change hebben zijn blauw gekleurd. Daarnaast zijn alle genen die geen statistisch significante p-waarde hebben en ook geen significante fold change hebben grijs gekleurd. De -log10(p-waarde) is op de y-as weergeven en op de x-as staat de log2(fold change).

Eigen analyse

Er zijn dus een paar genen te zien die wel interessant zouden zijn om te onderzoeken wat voor impact ze hebben op deze muizen. dit zijn bijvoorbeeld Mmp12, Hbb-bs of Fosb. Verder zijn de pax5 en ebf1 genen niet te zien in deze volcano plot.

Pathview:

Het probleem bij de eerste pathview is dat de expressies niveaus op een juiste manier hier in krijgen niet zo goed ging. Hierom heeft floris het gehele process nogmaals doorgelopen om te zien of hij de bugs kon fixen en het gehele process efficienter maken. Dit is geluk en er zijn een paar veranderingen. Floris heeft een pval_threshold van 0.05 hierdoor komen alleen statistisch significante resultaten in de deseq.degs.logfc

Van een aantal relevante pathways zijn afbeeldingen gemaakt, de Kegg pathway codes zijn afkomstig van de volgende website: https://www.genome.jp/kegg-bin/show_organism?menu_type=pathway_maps&org=mmu

De data die gebruikt wordt hier is de data van de deseq resultaten.

pval_threshold <- 0.05
data(gene.idtype.list); 
head(gene.idtype.list)
## [1] "SYMBOL"      "GENENAME"    "ENSEMBL"     "ENSEMBLPROT" "UNIGENE"    
## [6] "UNIPROT"
#hier word een subset gemaakt met de logfoldchange met 1 voorwaarde dat de padj waarde onder de pval_threshold ligt.
deseq.degs.logfc <- subset(res, padj < pval_threshold, select = log2FoldChange)

Het verschil met de pathway die ik gemaakt heb en die van floris is dus dat floris de logfc subset heeft en deze gebruikt als gene.data, in mijn pathview is de subset gebruikt

hier is nog een pathview voorbeeld van mij.

mmu04670_pathview <- pathview(gene.data = resultaat_deseq_datafr, pathway.id = “mmu04670”, species = “mmu”)

Als we deze vergelijken met die van floris is te zien dat ik de gehele subset mee geef en dat floris voor gene.data alleen een subset meegeeft met de logfoldchange waardoor de juiste expressie levels dus te zien zijn (-1 tot 1) ook is te zien dat er cpd.data= res toegevoegd is. De cpd staat voor compound en deze word terug gelinked aan de resultaten van de deseq.

Deze twee verandering hebben gezorgt dat de pathview er uit ziet zoals het moet en daarom zijn deze pathview resultaten leidend en gebruiken we deze voor de resultaten.

T cell receptor signaling pathway:

# 04660  T cell receptor signaling pathway

pathview(gene.data=deseq.degs.logfc,
        cpd.data=res, 
        ,gene.idtype = "SYMBOL",
        pathway.id="04660",
        species="mmu")
## 'select()' returned 1:1 mapping between keys and columns
## [1] "Note: 1 of 2119 unique input IDs unmapped."
## 'select()' returned 1:1 mapping between keys and columns
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04660.pathview.png
# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04660.pathview.png")
Figuur 27:  T cell receptor signaling pathway

Figuur 27: T cell receptor signaling pathway

T cell receptor signaling pathway: Er is overlap tussen de surface eiwitten van B- en T-lymfocyten. Dus ondanks dat de gesequencede cellen B-cellen zijn en geen T-cellen is het volgende pathway ook interesant. hier kunnen met over en onderexpressie van essentiele genen dus ook problemen ontstaan bij dingen zoals imuumsysteem reactie

# 04662  B cell receptor signaling pathway

pathview(gene.data=deseq.degs.logfc,
        cpd.data=res,
        ,gene.idtype = "SYMBOL",
        pathway.id="04662",
        species="mmu")
## 'select()' returned 1:1 mapping between keys and columns
## [1] "Note: 1 of 2119 unique input IDs unmapped."
## 'select()' returned 1:1 mapping between keys and columns
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04662.pathview.png
# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04662.pathview.png")
Figuur 28: B cell receptor signaling pathway

Figuur 28: B cell receptor signaling pathway

B cell receptor signaling pathway: De Genen die belangrijk zijn in ons originele onderzoek hebben afkomst van de B-cellen waardoor de B-cell receptoren pathway extreem relevant is. Wanneer er grote expressie verschillen zijn tussen de jonge en oude groep en er meer voorloper-b-cellen aanwezig zijn bij 1 van de groepen kan dit verklaard worden door expressie verschillen van de genen die verantwoordelijk zijn voor het maken van receptoreiwitten die als markers worden gebruikt. Hierdoor kunnen we ook zien dat wanneer er onder of overexpressie.

in onze pathway is te zien dat er een paar genen over expresses zijn en andere underexpressed. dit zou dus impact kunnen hebben op hoe sommige signalen worden verstuurd. Dit zou voor verschillende problemen kunnen zorgen. zoals problemen in imuum reactie.

# 04211

pathview(gene.data=deseq.degs.logfc,
        cpd.data=res,
        ,gene.idtype = "SYMBOL",
        pathway.id="04211",
        species="mmu")
## 'select()' returned 1:1 mapping between keys and columns
## [1] "Note: 1 of 2119 unique input IDs unmapped."
## 'select()' returned 1:1 mapping between keys and columns
## Warning: None of the genes or compounds mapped to the pathway!
## Argument gene.idtype or cpd.idtype may be wrong.
## Warning in colnames(plot.data)[c(1, 3, 9:ncs)] <- c("kegg.names", "all.mapped",
## : number of items to replace is not a multiple of replacement length
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04211.pathview.png
# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04211.pathview.png")
Figuur 29: Longevity regulating pathway

Figuur 29: Longevity regulating pathway

Analyse van floris

Longevity regulating pathway: De twee groepen die vergeleken worden zijn jonge en oude muizen, daarom is het longevity pathway dat enkele genen weergeeft die betrokken zijn bij het verouderingsproces ook belangrijk. Niet alle genen die betrokken zijn bij dit proces zijn hierop afgebeeld omdat het niet bekend is welke genen verantwoordelijk zijn. Foxo is downregulated in de oude muizengroep tegenover de jonge groep. Van dit gen is bekend dat het levensduur gerelateerd is. Bron: https://cellandbioscience.biomedcentral.com/articles/10.1186/s13578-021-00700-7

Analyse van mij

Problemen in deze pathway zullen dus impact hebben op de levensduur van de in de geval muis. Al is dit niet met zekerheid te zeggen omdat het nog niet bekend is welke genen specifiek belangrijk zijn voor levensduur.

# 04640

pathview(gene.data=deseq.degs.logfc,
        cpd.data=res,
        ,gene.idtype = "SYMBOL",
        pathway.id="04640",
        species="mmu")
## 'select()' returned 1:1 mapping between keys and columns
## [1] "Note: 1 of 2119 unique input IDs unmapped."
## 'select()' returned 1:1 mapping between keys and columns
## Info: Working in directory /Users/jarnoduiker/github_bioinf/Genomics_transcriptomics_analysis/analysis/Jarno
## Info: Writing image file mmu04640.pathview.png
# Om de afbeelding te weergeven in R:
knitr::include_graphics("mmu04640.pathview.png")
Figuur 30: Hematopoietic cell lineage

Figuur 30: Hematopoietic cell lineage

Hematopoietic cell lineage: Omdat B-cellen uit hematopoëtische cellen voortkwamen, en de voorloper B-cellen die hiervoor gebruikt zijn nog dichter in de buurt van deze hematopoëtische stamcellen zijn is het volgende pathway ook interessant. Deze pathway maakt ook verschillende dingen zoals de regulatory T-cell en als hier onderexpressie van in zou dit dus invloed hebben op regulatie in de T-cellen. Deze pathway heeft invloed op heel veel essentiele cellen maar omdat nog niet alles bekend is is het moeilijk om te zeggen wat deze over / onder expressies precies doen.

Resultaten & conclusie

Als we kijken naar dit onderzoek hebben we verschillende dingen bekeken, De pathways, De expressie van Ebf1 + Pax5 in jonge en oude muizen en andere genen die hoge expressie hadden. De terug koppeling naar ons originele onderzoek over de verschill van expressie tussen jong en oud van de dus essentiele Ebf1 en Pax5 genen is dat uit onze data is gebleken dat oude muizen lage expressie van deze genen hebben en de jonge muizen juist een hoge expressie, waardoor ik speculeer dat het onderzoek het gelijk had over de A B compartiment (In compartiment A vindt actieve transcriptie plaats en in B niet of in mindere mate.Van Ebf1 is wel bekend dat het met leeftijd van compartiment A naar B gaat, dit komt ook overeen met de verlaagde expressie in de oude muizen (Ebf1) ). Dit zijn compartimenten in de nucleus waar de chromatine heen verplaats waardoor het niet of minder word expressed. In het originele onderzoek was er weinig gespeculeerd over conformatieveranderingen van het DNA dat de onderzochte genen bevat.Wel was experimenteel vastgesteld dat Ebf1 van compartiment A naar B veranderd, op basis daarvan heb ik het bovenstaande vermoeden geconstrueerd.

Er is verder weinig nieuwe informatie gevonden vergeleken met het originele onderzoek ondanks het nieuwe referentie genoom. Veel andere analyses die voor het artikel uitgevoerd waren, waren niet van toepassing voor dit onderzoek omdat het analyses betrok die hier niet gedaan zijn zoals: HI-C, VDJ-seq en HiChIP.